]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/GammaConv/AliAnalysisTaskGammaConversion.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Changes from RR for PbPb, correction fromMW for direct photon
[u/mrichter/AliRoot.git] / PWG4 / GammaConv / AliAnalysisTaskGammaConversion.cxx
CommitLineData
d7d7e825 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: Ana Marin, Kathrin Koch, Kenneth Aamodt *
5 * Version 1.1 *
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////////////////////////////////////////////////
17//---------------------------------------------
18// Class used to do analysis on conversion pairs
19//---------------------------------------------
20////////////////////////////////////////////////
21
22// root
23#include <TChain.h>
24
25// analysis
26#include "AliAnalysisTaskGammaConversion.h"
27#include "AliStack.h"
28#include "AliLog.h"
29#include "AliESDtrackCuts.h"
30#include "TNtuple.h"
4a6157dc 31//#include "AliCFManager.h" // for CF
32//#include "AliCFContainer.h" // for CF
9c1cb6f7 33#include "AliESDInputHandler.h"
34#include "AliAnalysisManager.h"
6441e967 35#include "AliAODPWG4Particle.h"
36#include "AliAODPWG4ParticleCorrelation.h"
4a6157dc 37#include "AliGammaConversionAODObject.h"
6441e967 38#include "AliAODConversionParticle.h"
5e55d806 39#include "AliGammaConversionBGHandler.h"
6272370b 40#include "AliESDCaloCluster.h" // for combining PHOS and GammaConv
9c1cb6f7 41#include "AliKFVertex.h"
3c45d101 42#include "AliGenPythiaEventHeader.h"
43#include "AliGenDPMjetEventHeader.h"
44#include "AliGenEventHeader.h"
e60f3265 45#include <AliMCEventHandler.h>
77ac6f3e 46#include "TRandom3.h"
18a1e325 47#include "AliTriggerAnalysis.h"
b6dd6ad2 48#include "AliCentrality.h"
77ac6f3e 49
d707e3cf 50class AliESDTrackCuts;
4a6157dc 51class AliCFContainer;
52class AliCFManager;
d7d7e825 53class AliKFVertex;
54class AliAODHandler;
55class AliAODEvent;
56class ALiESDEvent;
57class AliMCEvent;
58class AliMCEventHandler;
59class AliESDInputHandler;
60class AliAnalysisManager;
61class Riostream;
62class TFile;
63class TInterpreter;
64class TSystem;
65class TROOT;
66
67ClassImp(AliAnalysisTaskGammaConversion)
68
69
70AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion():
71AliAnalysisTaskSE(),
72 fV0Reader(NULL),
73 fStack(NULL),
a0b94e5c 74 fMCTruth(NULL), // for CF
4a6157dc 75 fGCMCEvent(NULL), // for CF
d7d7e825 76 fESDEvent(NULL),
77 fOutputContainer(NULL),
a0b94e5c 78 fCFManager(0x0), // for CF
d7d7e825 79 fHistograms(NULL),
b5832f95 80 fTriggerCINT1B(kFALSE),
d7d7e825 81 fDoMCTruth(kFALSE),
82 fDoNeutralMeson(kFALSE),
6272370b 83 fDoOmegaMeson(kFALSE),
d7d7e825 84 fDoJet(kFALSE),
85 fDoChic(kFALSE),
9c1cb6f7 86 fRecalculateV0ForGamma(kFALSE),
6c84d371 87 fKFReconstructedGammasTClone(NULL),
6272370b 88 fKFReconstructedPi0sTClone(NULL),
9c1cb6f7 89 fKFRecalculatedGammasTClone(NULL),
6c84d371 90 fCurrentEventPosElectronTClone(NULL),
91 fCurrentEventNegElectronTClone(NULL),
92 fKFReconstructedGammasCutTClone(NULL),
93 fPreviousEventTLVNegElectronTClone(NULL),
94 fPreviousEventTLVPosElectronTClone(NULL),
d7d7e825 95 fElectronv1(),
96 fElectronv2(),
6272370b 97 fGammav1(),
98 fGammav2(),
9c1cb6f7 99 fElectronRecalculatedv1(),
100 fElectronRecalculatedv2(),
d7d7e825 101 fElectronMass(-1),
102 fGammaMass(-1),
103 fPi0Mass(-1),
104 fEtaMass(-1),
105 fGammaWidth(-1),
106 fPi0Width(-1),
107 fEtaWidth(-1),
108 fMinOpeningAngleGhostCut(0.),
9640a3d1 109 fEsdTrackCuts(NULL),
d7d7e825 110 fCalculateBackground(kFALSE),
111 fWriteNtuple(kFALSE),
112 fGammaNtuple(NULL),
113 fNeutralMesonNtuple(NULL),
114 fTotalNumberOfAddedNtupleEntries(0),
6c84d371 115 fChargedParticles(NULL),
d7d7e825 116 fChargedParticlesId(),
117 fGammaPtHighest(0.),
118 fMinPtForGammaJet(1.),
119 fMinIsoConeSize(0.2),
120 fMinPtIsoCone(0.7),
121 fMinPtGamChargedCorr(0.5),
122 fMinPtJetCone(0.5),
123 fLeadingChargedIndex(-1),
9640a3d1 124 fLowPtMapping(1.),
125 fHighPtMapping(3.),
1e7846f4 126 fDoCF(kFALSE),
04bf4381 127 fAODGamma(NULL),
128 fAODPi0(NULL),
129 fAODOmega(NULL),
037dc2db 130 fAODBranchName("GammaConv"),
6441e967 131 fOutputAODClassName("AliAODConversionParticle"),
63477d82 132 fKFCreateAOD(kTRUE),
d765d400 133 fKFForceAOD(kFALSE),
332f1f44 134 fKFDeltaAODFileName(""),
037dc2db 135 fDoNeutralMesonV0MCCheck(kFALSE),
5ce758b0 136 fUseTrackMultiplicityForBG(kTRUE),
137 fMoveParticleAccordingToVertex(kFALSE),
2e2da371 138 fApplyChi2Cut(kFALSE),
77ac6f3e 139 fNRandomEventsForBG(15),
140 fNDegreesPMBackground(15),
111d75df 141 fDoRotation(kTRUE),
142 fCheckBGProbability(kTRUE),
2f426a53 143 fKFReconstructedGammasV0Index(),
18a1e325 144 fRemovePileUp(kFALSE),
145 fSelectV0AND(kFALSE),
6746e1e1 146 fTriggerAnalysis(NULL),
147 fMultiplicity(0),
148 fUseMultiplicity(0),
64f4118c 149 fUseMultiplicityBin(0),
150 fUseCentrality(0),
151 fUseCentralityBin(0)
d7d7e825 152{
153 // Default constructor
5a34881d 154
155 /* Kenneth: the default constructor should not have any define input/output or the call to SetESDtrackCuts
d7d7e825 156 // Common I/O in slot 0
157 DefineInput (0, TChain::Class());
158 DefineOutput(0, TTree::Class());
159
160 // Your private output
161 DefineOutput(1, TList::Class());
a0b94e5c 162
d7d7e825 163 // Define standard ESD track cuts for Gamma-hadron correlation
164 SetESDtrackCuts();
5a34881d 165 */
d7d7e825 166}
167
168AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):
169 AliAnalysisTaskSE(name),
170 fV0Reader(NULL),
171 fStack(NULL),
a0b94e5c 172 fMCTruth(NULL), // for CF
4a6157dc 173 fGCMCEvent(NULL), // for CF
d7d7e825 174 fESDEvent(NULL),
175 fOutputContainer(0x0),
a0b94e5c 176 fCFManager(0x0), // for CF
d7d7e825 177 fHistograms(NULL),
b5832f95 178 fTriggerCINT1B(kFALSE),
d7d7e825 179 fDoMCTruth(kFALSE),
180 fDoNeutralMeson(kFALSE),
6272370b 181 fDoOmegaMeson(kFALSE),
d7d7e825 182 fDoJet(kFALSE),
183 fDoChic(kFALSE),
9c1cb6f7 184 fRecalculateV0ForGamma(kFALSE),
6c84d371 185 fKFReconstructedGammasTClone(NULL),
6272370b 186 fKFReconstructedPi0sTClone(NULL),
9c1cb6f7 187 fKFRecalculatedGammasTClone(NULL),
6c84d371 188 fCurrentEventPosElectronTClone(NULL),
189 fCurrentEventNegElectronTClone(NULL),
190 fKFReconstructedGammasCutTClone(NULL),
191 fPreviousEventTLVNegElectronTClone(NULL),
192 fPreviousEventTLVPosElectronTClone(NULL),
d7d7e825 193 fElectronv1(),
194 fElectronv2(),
6272370b 195 fGammav1(),
196 fGammav2(),
9c1cb6f7 197 fElectronRecalculatedv1(),
198 fElectronRecalculatedv2(),
d7d7e825 199 fElectronMass(-1),
200 fGammaMass(-1),
201 fPi0Mass(-1),
202 fEtaMass(-1),
203 fGammaWidth(-1),
204 fPi0Width(-1),
205 fEtaWidth(-1),
206 fMinOpeningAngleGhostCut(0.),
9640a3d1 207 fEsdTrackCuts(NULL),
d7d7e825 208 fCalculateBackground(kFALSE),
209 fWriteNtuple(kFALSE),
210 fGammaNtuple(NULL),
211 fNeutralMesonNtuple(NULL),
212 fTotalNumberOfAddedNtupleEntries(0),
6c84d371 213 fChargedParticles(NULL),
d7d7e825 214 fChargedParticlesId(),
215 fGammaPtHighest(0.),
216 fMinPtForGammaJet(1.),
217 fMinIsoConeSize(0.2),
218 fMinPtIsoCone(0.7),
219 fMinPtGamChargedCorr(0.5),
220 fMinPtJetCone(0.5),
221 fLeadingChargedIndex(-1),
9640a3d1 222 fLowPtMapping(1.),
223 fHighPtMapping(3.),
c59360eb 224 fDoCF(kFALSE),
225 fAODGamma(NULL),
226 fAODPi0(NULL),
227 fAODOmega(NULL),
228 fAODBranchName("GammaConv"),
6441e967 229 fOutputAODClassName("AliAODConversionParticle"),
63477d82 230 fKFCreateAOD(kTRUE),
c59360eb 231 fKFForceAOD(kFALSE),
232 fKFDeltaAODFileName(""),
233 fDoNeutralMesonV0MCCheck(kFALSE),
5ce758b0 234 fUseTrackMultiplicityForBG(kTRUE),
235 fMoveParticleAccordingToVertex(kFALSE),
2e2da371 236 fApplyChi2Cut(kFALSE),
77ac6f3e 237 fNRandomEventsForBG(15),
238 fNDegreesPMBackground(15),
111d75df 239 fDoRotation(kTRUE),
240 fCheckBGProbability(kTRUE),
2f426a53 241 fKFReconstructedGammasV0Index(),
18a1e325 242 fRemovePileUp(kFALSE),
243 fSelectV0AND(kFALSE),
6746e1e1 244 fTriggerAnalysis(NULL),
245 fMultiplicity(0),
246 fUseMultiplicity(0),
64f4118c 247 fUseMultiplicityBin(0),
248 fUseCentrality(0),
249 fUseCentralityBin(0)
d7d7e825 250{
251 // Common I/O in slot 0
252 DefineInput (0, TChain::Class());
253 DefineOutput(0, TTree::Class());
254
255 // Your private output
256 DefineOutput(1, TList::Class());
a0b94e5c 257 DefineOutput(2, AliCFContainer::Class()); // for CF
258
259
d7d7e825 260 // Define standard ESD track cuts for Gamma-hadron correlation
261 SetESDtrackCuts();
9c1cb6f7 262
d7d7e825 263}
264
265AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion()
266{
267 // Remove all pointers
268
269 if(fOutputContainer){
270 fOutputContainer->Clear() ;
271 delete fOutputContainer ;
272 }
273 if(fHistograms){
274 delete fHistograms;
275 }
276 if(fV0Reader){
277 delete fV0Reader;
278 }
a0b94e5c 279
280 // for CF
281 if(fCFManager){
282 delete fCFManager;
283 }
9640a3d1 284
285 if(fEsdTrackCuts){
286 delete fEsdTrackCuts;
287 }
288
04bf4381 289 //Delete AODs
290 if (fAODGamma) {
291 fAODGamma->Clear();
292 delete fAODGamma;
d7d7e825 293 }
04bf4381 294 fAODGamma = NULL;
295
296 if (fAODPi0) {
297 fAODPi0->Clear();
298 delete fAODPi0;
299 }
300 fAODPi0 = NULL;
d7d7e825 301
04bf4381 302 if (fAODOmega) {
303 fAODOmega->Clear();
304 delete fAODOmega;
305 }
306 fAODOmega = NULL;
307
18a1e325 308 if(fTriggerAnalysis) {
309 delete fTriggerAnalysis;
310 }
311
312
04bf4381 313}
d7d7e825 314
12464034 315
d7d7e825 316void AliAnalysisTaskGammaConversion::Init()
317{
318 // Initialization
4a6157dc 319 // AliLog::SetGlobalLogLevel(AliLog::kError);
d7d7e825 320}
321void AliAnalysisTaskGammaConversion::SetESDtrackCuts()
322{
323 // SetESDtrackCuts
9640a3d1 324 if (fEsdTrackCuts!=NULL){
325 delete fEsdTrackCuts;
326 }
d7d7e825 327 fEsdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
a0b94e5c 328 //standard cuts from:
329 //http://aliceinfo.cern.ch/alicvs/viewvc/PWG0/dNdEta/CreateCuts.C?revision=1.4&view=markup
ebcfaa7e 330
331 // Cuts used up to 3rd of March
332
333 // fEsdTrackCuts->SetMinNClustersTPC(50);
334// fEsdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
335// fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
336// fEsdTrackCuts->SetRequireITSRefit(kTRUE);
337// fEsdTrackCuts->SetMaxNsigmaToVertex(3);
338// fEsdTrackCuts->SetRequireSigmaToVertex(kTRUE);
339
340 //------- To be tested-----------
d707e3cf 341 // Cuts used up to 26th of Agost
342// Int_t minNClustersTPC = 70;
343// Double_t maxChi2PerClusterTPC = 4.0;
344// Double_t maxDCAtoVertexXY = 2.4; // cm
345// Double_t maxDCAtoVertexZ = 3.2; // cm
346// fEsdTrackCuts->SetRequireSigmaToVertex(kFALSE);
347// fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
348// fEsdTrackCuts->SetRequireITSRefit(kTRUE);
349// // fEsdTrackCuts->SetRequireTPCStandAlone(kTRUE);
350// fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
351// fEsdTrackCuts->SetMinNClustersTPC(minNClustersTPC);
352// fEsdTrackCuts->SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC);
353// fEsdTrackCuts->SetMaxDCAToVertexXY(maxDCAtoVertexXY);
354// fEsdTrackCuts->SetMaxDCAToVertexZ(maxDCAtoVertexZ);
355// fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
356// fEsdTrackCuts->SetEtaRange(-0.8, 0.8);
357// fEsdTrackCuts->SetPtRange(0.15);
358
359// fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst);
360
361
362// Using standard function for setting Cuts
363 Bool_t selectPrimaries=kTRUE;
406a77e0 364 fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
64f4118c 365 fEsdTrackCuts->SetMaxDCAToVertexZ(2);
d707e3cf 366 fEsdTrackCuts->SetEtaRange(-0.8, 0.8);
367 fEsdTrackCuts->SetPtRange(0.15);
ebcfaa7e 368
369 //----- From Jacek 10.03.03 ------------------/
370// minNClustersTPC = 70;
371// maxChi2PerClusterTPC = 4.0;
372// maxDCAtoVertexXY = 2.4; // cm
373// maxDCAtoVertexZ = 3.2; // cm
374
375// esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
376// esdTrackCuts->SetRequireTPCRefit(kFALSE);
377// esdTrackCuts->SetRequireTPCStandAlone(kTRUE);
378// esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
379// esdTrackCuts->SetMinNClustersTPC(minNClustersTPC);
380// esdTrackCuts->SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC);
381// esdTrackCuts->SetMaxDCAToVertexXY(maxDCAtoVertexXY);
382// esdTrackCuts->SetMaxDCAToVertexZ(maxDCAtoVertexZ);
383// esdTrackCuts->SetDCAToVertex2D(kTRUE);
384
385
386
d7d7e825 387 // fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
037dc2db 388 // fV0Reader->SetESDtrackCuts(fEsdTrackCuts);
d7d7e825 389}
390
c00009fb 391void AliAnalysisTaskGammaConversion::UserExec(Option_t */*option*/)
d7d7e825 392{
393 // Execute analysis for current event
6272370b 394
9c1cb6f7 395 // Load the esdpid from the esdhandler if exists (tender was applied) otherwise set the Bethe Bloch parameters
3c45d101 396 Int_t eventQuality=-1;
9c1cb6f7 397
398 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
399 AliESDInputHandler *esdHandler=0x0;
400 if ( (esdHandler=dynamic_cast<AliESDInputHandler*>(man->GetInputEventHandler())) && esdHandler->GetESDpid() ){
401 AliV0Reader::SetESDpid(esdHandler->GetESDpid());
402 } else {
403 //load esd pid bethe bloch parameters depending on the existance of the MC handler
404 // yes: MC parameters
405 // no: data parameters
406 if (!AliV0Reader::GetESDpid()){
407 if (fMCEvent ) {
408 AliV0Reader::InitESDpid();
409 } else {
410 AliV0Reader::InitESDpid(1);
411 }
412 }
413 }
414
d765d400 415 //Must set fForceAOD to true for the AOD to get filled. Should only be done when running independent chain / train.
416 if(fKFForceAOD) {
417 if (!AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()) AliFatal("Cannot run ESD filter without an output event handler");
418 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
419 }
48682642 420
421 if(fV0Reader == NULL){
422 // Write warning here cuts and so on are default if this ever happens
6272370b 423 }
48682642 424
e60f3265 425 if (fMCEvent ) {
426 // To avoid crashes due to unzip errors. Sometimes the trees are not there.
427
428 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
3c45d101 429 if (!mcHandler){
430 AliError("Could not retrive MC event handler!");
431 return;
432
433 eventQuality=0;
434 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
435 }
436 if (!mcHandler->InitOk() ){
437 eventQuality=0;
438 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
439 return;
440 }
441 if (!mcHandler->TreeK() ){
442 eventQuality=0;
443 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
444 return;
445 }
446 if (!mcHandler->TreeTR() ) {
447 eventQuality=0;
448 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
449 return;
450 }
e60f3265 451 }
452
48682642 453 fV0Reader->SetInputAndMCEvent(InputEvent(), MCEvent());
454
455 fV0Reader->Initialize();
456 fDoMCTruth = fV0Reader->GetDoMCTruth();
457
12464034 458 if(fAODGamma) fAODGamma->Delete();
459 if(fAODPi0) fAODPi0->Delete();
460 if(fAODOmega) fAODOmega->Delete();
a0b94e5c 461
6c84d371 462 if(fKFReconstructedGammasTClone == NULL){
463 fKFReconstructedGammasTClone = new TClonesArray("AliKFParticle",0);
464 }
465 if(fCurrentEventPosElectronTClone == NULL){
466 fCurrentEventPosElectronTClone = new TClonesArray("AliESDtrack",0);
467 }
468 if(fCurrentEventNegElectronTClone == NULL){
469 fCurrentEventNegElectronTClone = new TClonesArray("AliESDtrack",0);
470 }
471 if(fKFReconstructedGammasCutTClone == NULL){
472 fKFReconstructedGammasCutTClone = new TClonesArray("AliKFParticle",0);
473 }
474 if(fPreviousEventTLVNegElectronTClone == NULL){
475 fPreviousEventTLVNegElectronTClone = new TClonesArray("TLorentzVector",0);
476 }
477 if(fPreviousEventTLVPosElectronTClone == NULL){
478 fPreviousEventTLVPosElectronTClone = new TClonesArray("TLorentzVector",0);
479 }
480 if(fChargedParticles == NULL){
481 fChargedParticles = new TClonesArray("AliESDtrack",0);
482 }
6272370b 483
484 if(fKFReconstructedPi0sTClone == NULL){
485 fKFReconstructedPi0sTClone = new TClonesArray("AliKFParticle",0);
486 }
9c1cb6f7 487
18a1e325 488 if(fKFRecalculatedGammasTClone == NULL){
9c1cb6f7 489 fKFRecalculatedGammasTClone = new TClonesArray("AliKFParticle",0);
490 }
491
18a1e325 492 if(fTriggerAnalysis== NULL){
493 fTriggerAnalysis = new AliTriggerAnalysis;
494 }
495
6c84d371 496 //clear TClones
1e7846f4 497 fKFReconstructedGammasTClone->Delete();
498 fCurrentEventPosElectronTClone->Delete();
499 fCurrentEventNegElectronTClone->Delete();
500 fKFReconstructedGammasCutTClone->Delete();
501 fPreviousEventTLVNegElectronTClone->Delete();
502 fPreviousEventTLVPosElectronTClone->Delete();
6272370b 503 fKFReconstructedPi0sTClone->Delete();
9c1cb6f7 504 fKFRecalculatedGammasTClone->Delete();
5e55d806 505
d7d7e825 506 //clear vectors
d7d7e825 507 fElectronv1.clear();
508 fElectronv2.clear();
6272370b 509
510 fGammav1.clear();
511 fGammav2.clear();
512
9c1cb6f7 513 fElectronRecalculatedv1.clear();
514 fElectronRecalculatedv2.clear();
515
d7d7e825 516
1e7846f4 517 fChargedParticles->Delete();
5e55d806 518
d7d7e825 519 fChargedParticlesId.clear();
037dc2db 520
521 fKFReconstructedGammasV0Index.clear();
a0b94e5c 522
d7d7e825 523 //Clear the data in the v0Reader
5e55d806 524 // fV0Reader->UpdateEventByEventData();
b5832f95 525
526 //Take Only events with proper trigger
5e55d806 527 /*
b5832f95 528 if(fTriggerCINT1B){
529 if(!fV0Reader->GetESDEvent()->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) return;
530 }
5e55d806 531 */
22182c37 532 Bool_t v0A = fTriggerAnalysis->IsOfflineTriggerFired(fV0Reader->GetESDEvent(), AliTriggerAnalysis::kV0A);
533 Bool_t v0C = fTriggerAnalysis->IsOfflineTriggerFired(fV0Reader->GetESDEvent(), AliTriggerAnalysis::kV0C);
534 Bool_t v0AND = v0A && v0C;
535
536 if(fSelectV0AND && !v0AND){
537 eventQuality=5;
538 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
539 if(fDoMCTruth){
540 CheckMesonProcessTypeEventQuality(eventQuality);
541 }
542
543 return;
544 }
10e3319b 545
c8206114 546 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
547 // cout<< "Event not taken"<< endl;
3c45d101 548 eventQuality=1;
549 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
550 if(fDoMCTruth){
551 CheckMesonProcessTypeEventQuality(eventQuality);
552 }
c8206114 553 return; // aborts if the primary vertex does not have contributors.
554 }
555
ada7a878 556
c8206114 557
885114d1 558 if(!fV0Reader->CheckForPrimaryVertexZ() ){
3c45d101 559 eventQuality=2;
560 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
561 if(fDoMCTruth){
562 CheckMesonProcessTypeEventQuality(eventQuality);
563 }
885114d1 564 return;
565 }
2f426a53 566
ada7a878 567 if(fV0Reader->GetESDEvent()->GetPrimaryVertexTracks()->GetNContributors()>0) {
568 fHistograms->FillHistogram("ESD_GlobalPrimaryVtxZ",fV0Reader->GetESDEvent()->GetPrimaryVertex()->GetZ());
569 }else{
570 if(fV0Reader->GetESDEvent()->GetPrimaryVertexSPD()->GetNContributors()>0) {
571 fHistograms->FillHistogram("ESD_SPDPrimaryVtxZ",fV0Reader->GetESDEvent()->GetPrimaryVertex()->GetZ());
572 }
573 }
574
2f426a53 575 if(fRemovePileUp && fV0Reader->GetESDEvent()->IsPileupFromSPD()) {
0a3d9f49 576 eventQuality=4;
45441b78 577 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
2f426a53 578 return;
579 }
18a1e325 580
6746e1e1 581 fMultiplicity = fEsdTrackCuts->CountAcceptedTracks(fV0Reader->GetESDEvent());
582
583 if( CalculateMultiplicityBin() != fUseMultiplicityBin){
584 eventQuality=6;
585 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
586 return;
587 }
18a1e325 588
64f4118c 589 if(fV0Reader->GetIsHeavyIon()){
590 if(fUseCentrality>0){
b6dd6ad2 591 AliCentrality *esdCentrality = fV0Reader->GetESDEvent()->GetCentrality();
64f4118c 592 Int_t centralityC = -1;
593
594 if(fUseCentrality==1){
595 centralityC = esdCentrality->GetCentralityClass10("V0M");
596 if( centralityC != fUseCentralityBin ){
597 eventQuality=7;
598 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
599 return;
600 }
601 }
602
603 if(fUseCentrality==2){
604 centralityC = esdCentrality->GetCentralityClass10("CL1");
605 if( centralityC != fUseCentralityBin ){
606 eventQuality=7;
607 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
608 return;
609 }
610 }
14d7f2f9 611
612////////////////////////////////////// RRnew start /////////////////////////////////////////////////////
613 if(fUseCentrality==3){
614 centralityC = esdCentrality->GetCentralityClass10("V0M");
615 if( (fUseCentralityBin == 0) && (centralityC!=0) ){ // 0-10%
616 eventQuality=7;
617 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
618 return;
619 }
620 if( (fUseCentralityBin == 1) && (centralityC!=1) ){ // 10-20%
621 eventQuality=7;
622 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
623 return;
624 }
625 if( (fUseCentralityBin == 2) && (centralityC!=2) && (centralityC!=3) ){ // 20-40%
626 eventQuality=7;
627 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
628 return;
629 }
630 if( (fUseCentralityBin == 4) && (centralityC!=4) && (centralityC!=5) ){ // 40-60%
631 eventQuality=7;
632 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
633 return;
634 }
635 if( (fUseCentralityBin == 6) && (centralityC!=6) && (centralityC!=7) && (centralityC!=8) ){ // 60-90%
636 eventQuality=7;
637 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
638 return;
639 }
640 }
641
642 if(fUseCentrality==4){
643 centralityC = esdCentrality->GetCentralityClass10("CL1");
644 if( (fUseCentralityBin == 0) && (centralityC!=0) ){ // 0-10%
645 eventQuality=7;
646 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
647 return;
648 }
649 if( (fUseCentralityBin == 1) && (centralityC!=1) ){ // 10-20%
650 eventQuality=7;
651 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
652 return;
653 }
654 if( (fUseCentralityBin == 2) && (centralityC!=2) && (centralityC!=3) ){ // 20-40%
655 eventQuality=7;
656 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
657 return;
658 }
659 if( (fUseCentralityBin == 4) && (centralityC!=4) && (centralityC!=5) ){ // 40-60%
660 eventQuality=7;
661 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
662 return;
663 }
664 if( (fUseCentralityBin == 6) && (centralityC!=6) && (centralityC!=7) && (centralityC!=8) ){ // 60-90%
665 eventQuality=7;
666 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
667 return;
668 }
669 }
670////////////////////////////////////// RRnew end ///////////////////////////////////////////////////////
671
64f4118c 672 }
673 }
3c45d101 674 eventQuality=3;
675 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
885114d1 676
6746e1e1 677
678
679 fHistograms->FillHistogram("ESD_NumberOfGoodESDTracks",fMultiplicity);
680 if (fV0Reader->GetNumberOfContributorsVtx()>=1){
681 fHistograms->FillHistogram("ESD_NumberOfGoodESDTracksVtx",fMultiplicity);
682 }
683
684
685
d7d7e825 686 // Process the MC information
687 if(fDoMCTruth){
688 ProcessMCData();
689 }
a0b94e5c 690
d7d7e825 691 //Process the v0 information with no cuts
692 ProcessV0sNoCut();
cb90a330 693
d7d7e825 694 // Process the v0 information
695 ProcessV0s();
cb90a330 696
5e55d806 697
d7d7e825 698 //Fill Gamma AOD
63477d82 699 if(fKFCreateAOD) {
700 FillAODWithConversionGammas() ;
701 }
702
703
d7d7e825 704 // Process reconstructed gammas
705 if(fDoNeutralMeson == kTRUE){
a0b94e5c 706 ProcessGammasForNeutralMesonAnalysis();
77880bd8 707
48682642 708 }
709
77880bd8 710 if(fDoMCTruth == kTRUE){
711 CheckV0Efficiency();
712 }
d7d7e825 713 //Process reconstructed gammas electrons for Chi_c Analysis
6c84d371 714 if(fDoChic == kTRUE){
d7d7e825 715 ProcessGammaElectronsForChicAnalysis();
716 }
717 // Process reconstructed gammas for gamma Jet/hadron correlations
718 if(fDoJet == kTRUE){
719 ProcessGammasForGammaJetAnalysis();
720 }
a0b94e5c 721
5e55d806 722 //calculate background if flag is set
723 if(fCalculateBackground){
724 CalculateBackground();
725 }
48682642 726
727 if(fDoNeutralMeson == kTRUE){
728// ProcessConvPHOSGammasForNeutralMesonAnalysis();
729 if(fDoOmegaMeson == kTRUE){
730 ProcessGammasForOmegaMesonAnalysis();
731 }
732 }
733
5e55d806 734 //Clear the data in the v0Reader
735 fV0Reader->UpdateEventByEventData();
9c1cb6f7 736 if(fRecalculateV0ForGamma==kTRUE){
737 RecalculateV0ForGamma();
738 }
d7d7e825 739 PostData(1, fOutputContainer);
a0b94e5c 740 PostData(2, fCFManager->GetParticleContainer()); // for CF
d7d7e825 741
742}
743
48682642 744// void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t *option){
745// // see header file for documentation
746// // printf(" ConnectInputData %s\n", GetName());
8a685cf3 747
48682642 748// AliAnalysisTaskSE::ConnectInputData(option);
8a685cf3 749
48682642 750// if(fV0Reader == NULL){
751// // Write warning here cuts and so on are default if this ever happens
752// }
753// fV0Reader->Initialize();
754// fDoMCTruth = fV0Reader->GetDoMCTruth();
755// }
d7d7e825 756
3c45d101 757void AliAnalysisTaskGammaConversion::CheckMesonProcessTypeEventQuality(Int_t evtQ){
77ac6f3e 758 // Check meson process type event quality
3c45d101 759 fStack= MCEvent()->Stack();
760 fGCMCEvent=MCEvent();
761
762 for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {
763 TParticle* particle = (TParticle *)fStack->Particle(iTracks);
764 if (!particle) {
765 //print warning here
766 continue;
767 }
768 if(particle->GetPdgCode()!=111){ //Pi0
769 continue;
770 }
62ae0124 771
772 Double_t rapidity;
773 if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){
774 rapidity=0;
775 }
776 else{
777 rapidity = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())));
778 }
779
780 if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ) continue;
f3026589 781 if(fV0Reader->GetIsHeavyIon()) continue;
62ae0124 782
3c45d101 783 if(evtQ==1){
784 switch(GetProcessType(fGCMCEvent)){
785 case kProcSD:
786 fHistograms->FillHistogram("MC_SD_EvtQ1_Pi0_Pt", particle->Pt());
787 break;
788 case kProcDD:
789 fHistograms->FillHistogram("MC_DD_EvtQ1_Pi0_Pt", particle->Pt());
790 break;
791 case kProcND:
792 fHistograms->FillHistogram("MC_ND_EvtQ1_Pi0_Pt", particle->Pt());
793 break;
794 default:
795 AliError("Unknown Process");
796 }
797 }
798 if(evtQ==2){
799 switch(GetProcessType(fGCMCEvent)){
800 case kProcSD:
801 fHistograms->FillHistogram("MC_SD_EvtQ2_Pi0_Pt", particle->Pt());
802 break;
803 case kProcDD:
804 fHistograms->FillHistogram("MC_DD_EvtQ2_Pi0_Pt", particle->Pt());
805 break;
806 case kProcND:
807 fHistograms->FillHistogram("MC_ND_EvtQ2_Pi0_Pt", particle->Pt());
808 break;
809 default:
810 AliError("Unknown Process");
811 }
812 }
d7d7e825 813
62ae0124 814 if(evtQ==4){
815 switch(GetProcessType(fGCMCEvent)){
816 case kProcSD:
817 fHistograms->FillHistogram("MC_SD_EvtQ4_Pi0_Pt", particle->Pt());
818 break;
819 case kProcDD:
820 fHistograms->FillHistogram("MC_DD_EvtQ4_Pi0_Pt", particle->Pt());
821 break;
822 case kProcND:
823 fHistograms->FillHistogram("MC_ND_EvtQ4_Pi0_Pt", particle->Pt());
824 break;
825 default:
826 AliError("Unknown Process");
827 }
828 }
829
830 if(evtQ==5){
831 switch(GetProcessType(fGCMCEvent)){
832 case kProcSD:
833 fHistograms->FillHistogram("MC_SD_EvtQ5_Pi0_Pt", particle->Pt());
834 break;
835 case kProcDD:
836 fHistograms->FillHistogram("MC_DD_EvtQ5_Pi0_Pt", particle->Pt());
837 break;
838 case kProcND:
839 fHistograms->FillHistogram("MC_ND_EvtQ5_Pi0_Pt", particle->Pt());
840 break;
841 default:
842 AliError("Unknown Process");
843 }
844 }
d7d7e825 845
3c45d101 846 }
847
848}
849
d7d7e825 850void AliAnalysisTaskGammaConversion::ProcessMCData(){
851 // see header file for documentation
48682642 852 //InputEvent(), MCEvent());
853 /* TestAnaMarin
d7d7e825 854 fStack = fV0Reader->GetMCStack();
a0b94e5c 855 fMCTruth = fV0Reader->GetMCTruth(); // for CF
4a6157dc 856 fGCMCEvent = fV0Reader->GetMCEvent(); // for CF
48682642 857 */
858 fStack= MCEvent()->Stack();
859 fGCMCEvent=MCEvent();
a0b94e5c 860
861 // for CF
1e7846f4 862 Double_t containerInput[3];
863 if(fDoCF){
864 if(!fGCMCEvent) cout << "NO MC INFO FOUND" << endl;
865 fCFManager->SetEventInfo(fGCMCEvent);
866 }
a0b94e5c 867 // end for CF
868
869
d7d7e825 870 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
871 return; // aborts if the primary vertex does not have contributors.
872 }
2f426a53 873 for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {
874 // for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {
d7d7e825 875 TParticle* particle = (TParticle *)fStack->Particle(iTracks);
bd6d9fa3 876
111d75df 877
bd6d9fa3 878
d7d7e825 879 if (!particle) {
880 //print warning here
881 continue;
882 }
a0b94e5c 883
d7d7e825 884 ///////////////////////Begin Chic Analysis/////////////////////////////
91b98865 885 if(fDoChic) {
886 if(particle->GetPdgCode() == 443){//Is JPsi
887 if(particle->GetNDaughters()==2){
888 if(TMath::Abs(fStack->Particle(particle->GetFirstDaughter())->GetPdgCode()) == 11 &&
889 TMath::Abs(fStack->Particle(particle->GetLastDaughter())->GetPdgCode()) == 11){
890
891 TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
892 TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
893 if(TMath::Abs(daug0->Eta()) < 0.9 && TMath::Abs(daug1->Eta()) < 0.9)
894 fHistograms->FillTable("Table_Electrons",3);//e+ e- from J/Psi inside acceptance
a0b94e5c 895
91b98865 896 if( TMath::Abs(daug0->Eta()) < 0.9){
897 if(daug0->GetPdgCode() == -11)
898 fHistograms->FillTable("Table_Electrons",1);//e+ from J/Psi inside acceptance
899 else
900 fHistograms->FillTable("Table_Electrons",2);//e- from J/Psi inside acceptance
a0b94e5c 901
91b98865 902 }
903 if(TMath::Abs(daug1->Eta()) < 0.9){
904 if(daug1->GetPdgCode() == -11)
905 fHistograms->FillTable("Table_Electrons",1);//e+ from J/Psi inside acceptance
906 else
907 fHistograms->FillTable("Table_Electrons",2);//e- from J/Psi inside acceptance
908 }
d7d7e825 909 }
910 }
911 }
91b98865 912 // const int CHI_C0 = 10441;
913 // const int CHI_C1 = 20443;
914 // const int CHI_C2 = 445
915 if(particle->GetPdgCode() == 22){//gamma from JPsi
916 if(particle->GetMother(0) > -1){
917 if(fStack->Particle(particle->GetMother(0))->GetPdgCode() == 10441 ||
918 fStack->Particle(particle->GetMother(0))->GetPdgCode() == 20443 ||
919 fStack->Particle(particle->GetMother(0))->GetPdgCode() == 445){
920 if(TMath::Abs(particle->Eta()) < 1.2)
921 fHistograms->FillTable("Table_Electrons",17);// gamma from chic inside accptance
922 }
d7d7e825 923 }
924 }
91b98865 925 if(particle->GetPdgCode() == 10441 || particle->GetPdgCode() == 20443 || particle->GetPdgCode() == 445){
926 if( particle->GetNDaughters() == 2){
927 TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
928 TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
a0b94e5c 929
91b98865 930 if( (daug0->GetPdgCode() == 443 || daug0->GetPdgCode() == 22) && (daug1->GetPdgCode() == 443 || daug1->GetPdgCode() == 22) ){
931 if( daug0->GetPdgCode() == 443){
932 TParticle* daugE0 = fStack->Particle(daug0->GetFirstDaughter());
933 TParticle* daugE1 = fStack->Particle(daug0->GetLastDaughter());
934 if( TMath::Abs(daug1->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
935 fHistograms->FillTable("Table_Electrons",18);
a0b94e5c 936
91b98865 937 }//if
938 else if (daug1->GetPdgCode() == 443){
939 TParticle* daugE0 = fStack->Particle(daug1->GetFirstDaughter());
940 TParticle* daugE1 = fStack->Particle(daug1->GetLastDaughter());
941 if( TMath::Abs(daug0->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
942 fHistograms->FillTable("Table_Electrons",18);
943 }//else if
944 }//gamma o Jpsi
945 }//GetNDaughters
946 }
d7d7e825 947 }
a0b94e5c 948
d7d7e825 949 /////////////////////End Chic Analysis////////////////////////////
a0b94e5c 950
951
dc2883e4 952 // if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() ) continue;
a0b94e5c 953
d7d7e825 954 if(particle->R()>fV0Reader->GetMaxRCut()) continue; // cuts on distance from collision point
955
956 Double_t tmpPhi=particle->Phi();
957
958 if(particle->Phi()> TMath::Pi()){
959 tmpPhi = particle->Phi()-(2*TMath::Pi());
960 }
961
962 Double_t rapidity;
963 if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){
964 rapidity=0;
965 }
966 else{
967 rapidity = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())));
968 }
969
111d75df 970
971
972 if(iTracks<=fStack->GetNprimary() ){
973 if ( particle->GetPdgCode()== -211 || particle->GetPdgCode()== 211 ||
974 particle->GetPdgCode()== 2212 || particle->GetPdgCode()==-2212 ||
975 particle->GetPdgCode()== 321 || particle->GetPdgCode()==-321 ){
dc2883e4 976 if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() ) continue;
3c45d101 977 fHistograms->FillHistogram("MC_PhysicalPrimaryCharged_Pt", particle->Pt());
111d75df 978 }
bd6d9fa3 979 }
111d75df 980
bd6d9fa3 981
982
d7d7e825 983 //process the gammas
984 if (particle->GetPdgCode() == 22){
dc2883e4 985 if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() ) continue;
a607c218 986
d7d7e825 987 if(particle->GetMother(0) >-1 && fStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){
988 continue; // no photon as mothers!
989 }
a607c218 990
d7d7e825 991 if(particle->GetMother(0) >= fStack->GetNprimary()){
992 continue; // the gamma has a mother, and it is not a primary particle
993 }
f8017b04 994
995 if(particle->GetMother(0) >-1){
996 fHistograms->FillHistogram("MC_DecayAllGamma_Pt", particle->Pt()); // All
997 switch(fStack->Particle(particle->GetMother(0))->GetPdgCode()){
998 case 111: // Pi0
999 fHistograms->FillHistogram("MC_DecayPi0Gamma_Pt", particle->Pt());
1000 break;
1001 case 113: // Rho0
1002 fHistograms->FillHistogram("MC_DecayRho0Gamma_Pt", particle->Pt());
1003 break;
1004 case 221: // Eta
1005 fHistograms->FillHistogram("MC_DecayEtaGamma_Pt", particle->Pt());
1006 break;
1007 case 223: // Omega
1008 fHistograms->FillHistogram("MC_DecayOmegaGamma_Pt", particle->Pt());
1009 break;
1010 case 310: // K_s0
1011 fHistograms->FillHistogram("MC_DecayK0sGamma_Pt", particle->Pt());
1012 break;
1013 case 331: // Eta'
1014 fHistograms->FillHistogram("MC_DecayEtapGamma_Pt", particle->Pt());
1015 break;
86e91247 1016 case 333: // Phi
1017 fHistograms->FillHistogram("MC_DecayPhiGamma_Pt", particle->Pt());
1018 break;
f8017b04 1019 }
1020 }
1021
d7d7e825 1022 fHistograms->FillHistogram("MC_allGamma_Energy", particle->Energy());
1023 fHistograms->FillHistogram("MC_allGamma_Pt", particle->Pt());
1024 fHistograms->FillHistogram("MC_allGamma_Eta", particle->Eta());
1025 fHistograms->FillHistogram("MC_allGamma_Phi", tmpPhi);
1026 fHistograms->FillHistogram("MC_allGamma_Rapid", rapidity);
1027
a0b94e5c 1028 // for CF
1e7846f4 1029 if(fDoCF){
1030 containerInput[0] = particle->Pt();
1031 containerInput[1] = particle->Eta();
1032 if(particle->GetMother(0) >=0){
1033 containerInput[2] = fStack->Particle(particle->GetMother(0))->GetMass();
1034 }
1035 else{
1036 containerInput[2]=-1;
1037 }
a607c218 1038
1e7846f4 1039 fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated); // generated gamma
1040 }
86e91247 1041 if(particle->GetMother(0) < 0 || //Phojet p+p -> Direct Photons have no mother
1042 ((particle->GetMother(0) > -1) &&
14d7f2f9 1043 ((TMath::Abs(fStack->Particle(particle->GetMother(0))->GetPdgCode()) < 10)||
1044 (TMath::Abs(fStack->Particle(particle->GetMother(0))->GetPdgCode()) ==21) )) //Pythia p+p -> Direct Photons have quarksor gluons as mother
86e91247 1045 ){ // direct gamma
d7d7e825 1046 fHistograms->FillHistogram("MC_allDirectGamma_Energy",particle->Energy());
1047 fHistograms->FillHistogram("MC_allDirectGamma_Pt", particle->Pt());
1048 fHistograms->FillHistogram("MC_allDirectGamma_Eta", particle->Eta());
1049 fHistograms->FillHistogram("MC_allDirectGamma_Phi", tmpPhi);
1050 fHistograms->FillHistogram("MC_allDirectGamma_Rapid", rapidity);
1051 }
1052
d7d7e825 1053 // looking for conversion (electron + positron from pairbuilding (= 5) )
1054 TParticle* ePos = NULL;
1055 TParticle* eNeg = NULL;
1056
1057 if(particle->GetNDaughters() >= 2){
1058 for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
1059 TParticle *tmpDaughter = fStack->Particle(daughterIndex);
1060 if(tmpDaughter->GetUniqueID() == 5){
1061 if(tmpDaughter->GetPdgCode() == 11){
1062 eNeg = tmpDaughter;
1063 }
1064 else if(tmpDaughter->GetPdgCode() == -11){
1065 ePos = tmpDaughter;
1066 }
1067 }
1068 }
1069 }
1070
1071
1072 if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
1073 continue;
1074 }
1075
1076
1077 Double_t ePosPhi = ePos->Phi();
1078 if(ePos->Phi()> TMath::Pi()) ePosPhi = ePos->Phi()-(2*TMath::Pi());
1079
1080 Double_t eNegPhi = eNeg->Phi();
1081 if(eNeg->Phi()> TMath::Pi()) eNegPhi = eNeg->Phi()-(2*TMath::Pi());
1082
1083
1084 if(ePos->Pt()<fV0Reader->GetPtCut() || eNeg->Pt()<fV0Reader->GetPtCut()){
1085 continue; // no reconstruction below the Pt cut
1086 }
a0b94e5c 1087
d7d7e825 1088 if(TMath::Abs(ePos->Eta())> fV0Reader->GetEtaCut() || TMath::Abs(eNeg->Eta())> fV0Reader->GetEtaCut()){
1089 continue;
1090 }
a0b94e5c 1091
d7d7e825 1092 if(ePos->R()>fV0Reader->GetMaxRCut()){
1093 continue; // cuts on distance from collision point
1094 }
a0b94e5c 1095
1096 if(TMath::Abs(ePos->Vz()) > fV0Reader->GetMaxZCut()){
d7d7e825 1097 continue; // outside material
1098 }
1099
a0b94e5c 1100
d7d7e825 1101 if((TMath::Abs(ePos->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() > ePos->R()){
1102 continue; // line cut to exclude regions where we do not reconstruct
1103 }
1104
a0b94e5c 1105
1106 // for CF
1e7846f4 1107 if(fDoCF){
1108 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructable); // reconstructable gamma
1109 }
d7d7e825 1110 fHistograms->FillHistogram("MC_ConvGamma_Energy", particle->Energy());
1111 fHistograms->FillHistogram("MC_ConvGamma_Pt", particle->Pt());
1112 fHistograms->FillHistogram("MC_ConvGamma_Eta", particle->Eta());
1113 fHistograms->FillHistogram("MC_ConvGamma_Phi", tmpPhi);
1114 fHistograms->FillHistogram("MC_ConvGamma_Rapid", rapidity);
1115 fHistograms->FillHistogram("MC_ConvGamma_Pt_Eta", particle->Pt(),particle->Eta());
1116
1117 fHistograms->FillHistogram("MC_E_Energy", eNeg->Energy());
1118 fHistograms->FillHistogram("MC_E_Pt", eNeg->Pt());
1119 fHistograms->FillHistogram("MC_E_Eta", eNeg->Eta());
1120 fHistograms->FillHistogram("MC_E_Phi", eNegPhi);
1121
1122 fHistograms->FillHistogram("MC_P_Energy", ePos->Energy());
1123 fHistograms->FillHistogram("MC_P_Pt", ePos->Pt());
1124 fHistograms->FillHistogram("MC_P_Eta", ePos->Eta());
1125 fHistograms->FillHistogram("MC_P_Phi", ePosPhi);
1126
1127
d7d7e825 1128 // begin Mapping
1129 Int_t rBin = fHistograms->GetRBin(ePos->R());
9640a3d1 1130 Int_t zBin = fHistograms->GetZBin(ePos->Vz());
d7d7e825 1131 Int_t phiBin = fHistograms->GetPhiBin(particle->Phi());
9c1cb6f7 1132 Double_t rFMD=30;
64f4118c 1133 Double_t rITSTPCMin=50;
1134 Double_t rITSTPCMax=80;
9640a3d1 1135
1136 TVector3 vtxPos(ePos->Vx(),ePos->Vy(),ePos->Vz());
1137
d7d7e825 1138 TString nameMCMappingPhiR="";
e158cbc3 1139 nameMCMappingPhiR.Form("MC_Conversion_Mapping_Phi%02d_R%02d",phiBin,rBin);
ebcfaa7e 1140 // fHistograms->FillHistogram(nameMCMappingPhiR, ePos->Vz(), particle->Eta());
d7d7e825 1141
1142 TString nameMCMappingPhi="";
e158cbc3 1143 nameMCMappingPhi.Form("MC_Conversion_Mapping_Phi%02d",phiBin);
9640a3d1 1144 // fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());
ebcfaa7e 1145 //fHistograms->FillHistogram(nameMCMappingPhi, ePos->Vz(), particle->Eta());
d7d7e825 1146
1147 TString nameMCMappingR="";
e158cbc3 1148 nameMCMappingR.Form("MC_Conversion_Mapping_R%02d",rBin);
9640a3d1 1149 // fHistograms->FillHistogram(nameMCMappingR, particle->Eta());
ebcfaa7e 1150 //fHistograms->FillHistogram(nameMCMappingR,ePos->Vz(), particle->Eta());
d7d7e825 1151
1152 TString nameMCMappingPhiInR="";
e158cbc3 1153 nameMCMappingPhiInR.Form("MC_Conversion_Mapping_Phi_in_R_%02d",rBin);
9640a3d1 1154 // fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
1155 fHistograms->FillHistogram(nameMCMappingPhiInR, vtxPos.Phi());
1156
1157 TString nameMCMappingZInR="";
1158 nameMCMappingZInR.Form("MC_Conversion_Mapping_Z_in_R_%02d",rBin);
1159 fHistograms->FillHistogram(nameMCMappingZInR,ePos->Vz() );
1160
1161
1162 TString nameMCMappingPhiInZ="";
1163 nameMCMappingPhiInZ.Form("MC_Conversion_Mapping_Phi_in_Z_%02d",zBin);
1164 // fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
1165 fHistograms->FillHistogram(nameMCMappingPhiInZ, vtxPos.Phi());
1166
9c1cb6f7 1167
1168 if(ePos->R()<rFMD){
1169 TString nameMCMappingFMDPhiInZ="";
1170 nameMCMappingFMDPhiInZ.Form("MC_Conversion_Mapping_FMD_Phi_in_Z_%02d",zBin);
1171 fHistograms->FillHistogram(nameMCMappingFMDPhiInZ, vtxPos.Phi());
1172 }
1173
64f4118c 1174 if(ePos->R()>rITSTPCMin && ePos->R()<rITSTPCMax){
1175 TString nameMCMappingITSTPCPhiInZ="";
1176 nameMCMappingITSTPCPhiInZ.Form("MC_Conversion_Mapping_ITSTPC_Phi_in_Z_%02d",zBin);
1177 fHistograms->FillHistogram(nameMCMappingITSTPCPhiInZ, vtxPos.Phi());
1178 }
1179
9640a3d1 1180 TString nameMCMappingRInZ="";
1181 nameMCMappingRInZ.Form("MC_Conversion_Mapping_R_in_Z_%02d",zBin);
1182 fHistograms->FillHistogram(nameMCMappingRInZ,ePos->R() );
1183
1184 if(particle->Pt() > fLowPtMapping && particle->Pt()< fHighPtMapping){
1185 TString nameMCMappingMidPtPhiInR="";
1186 nameMCMappingMidPtPhiInR.Form("MC_Conversion_Mapping_MidPt_Phi_in_R_%02d",rBin);
1187 fHistograms->FillHistogram(nameMCMappingMidPtPhiInR, vtxPos.Phi());
1188
1189 TString nameMCMappingMidPtZInR="";
1190 nameMCMappingMidPtZInR.Form("MC_Conversion_Mapping_MidPt_Z_in_R_%02d",rBin);
1191 fHistograms->FillHistogram(nameMCMappingMidPtZInR,ePos->Vz() );
1192
1193
1194 TString nameMCMappingMidPtPhiInZ="";
1195 nameMCMappingMidPtPhiInZ.Form("MC_Conversion_Mapping_MidPt_Phi_in_Z_%02d",zBin);
1196 fHistograms->FillHistogram(nameMCMappingMidPtPhiInZ, vtxPos.Phi());
9c1cb6f7 1197
1198
1199 if(ePos->R()<rFMD){
1200 TString nameMCMappingMidPtFMDPhiInZ="";
1201 nameMCMappingMidPtFMDPhiInZ.Form("MC_Conversion_Mapping_MidPt_FMD_Phi_in_Z_%02d",zBin);
1202 fHistograms->FillHistogram(nameMCMappingMidPtFMDPhiInZ, vtxPos.Phi());
1203 }
1204
9640a3d1 1205
1206 TString nameMCMappingMidPtRInZ="";
1207 nameMCMappingMidPtRInZ.Form("MC_Conversion_Mapping_MidPt_R_in_Z_%02d",zBin);
1208 fHistograms->FillHistogram(nameMCMappingMidPtRInZ,ePos->R() );
1209
1210 }
1211
d7d7e825 1212 //end mapping
1213
1214 fHistograms->FillHistogram("MC_Conversion_R",ePos->R());
1215 fHistograms->FillHistogram("MC_Conversion_ZR",ePos->Vz(),ePos->R());
1216 fHistograms->FillHistogram("MC_Conversion_XY",ePos->Vx(),ePos->Vy());
1217 fHistograms->FillHistogram("MC_Conversion_OpeningAngle",GetMCOpeningAngle(ePos, eNeg));
9640a3d1 1218 fHistograms->FillHistogram("MC_ConvGamma_E_AsymmetryP",particle->P(),eNeg->P()/particle->P());
1219 fHistograms->FillHistogram("MC_ConvGamma_P_AsymmetryP",particle->P(),ePos->P()/particle->P());
1220
1221
d7d7e825 1222 if(particle->GetMother(0) < 0){ // no mother = direct gamma, still inside converted
1223 fHistograms->FillHistogram("MC_ConvDirectGamma_Energy",particle->Energy());
1224 fHistograms->FillHistogram("MC_ConvDirectGamma_Pt", particle->Pt());
1225 fHistograms->FillHistogram("MC_ConvDirectGamma_Eta", particle->Eta());
1226 fHistograms->FillHistogram("MC_ConvDirectGamma_Phi", tmpPhi);
1227 fHistograms->FillHistogram("MC_ConvDirectGamma_Rapid", rapidity);
1228
1229 } // end direct gamma
1230 else{ // mother exits
6c84d371 1231 /* if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0
a0b94e5c 1232 fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S
1233 fStack->Particle(particle->GetMother(0))->GetPdgCode()==445 //chic2
1234 ){
1235 fMCGammaChic.push_back(particle);
1236 }
6c84d371 1237 */
d7d7e825 1238 } // end if mother exits
1239 } // end if particle is a photon
1240
1241
1242
1243 // process motherparticles (2 gammas as daughters)
1244 // the motherparticle had already to pass the R and the eta cut, but no line cut.
1245 // the line cut is just valid for the conversions!
b8faf1bc 1246
14d7f2f9 1247 // RR primary Pi0 debug ////////////////////////////////////////////////////////////////////////////////////////////
b8faf1bc 1248 if (particle->GetPdgCode()==111){
1249 if( TMath::Abs(rapidity) < fV0Reader->GetRapidityMesonCut() ){
1250 fHistograms->FillHistogram("MC_Pi0_Pt_vs_Rapid_allDaughters", particle->Pt(),rapidity);
1251 }
1252 }
14d7f2f9 1253 // end RR primary Pi0 debug ////////////////////////////////////////////////////////////////////////////////////////
d7d7e825 1254
1255 if(particle->GetNDaughters() == 2){
1256
1257 TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());
1258 TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());
1259
1260 if(daughter0->GetPdgCode() != 22 || daughter1->GetPdgCode() != 22) continue; //check for gamma gamma daughters
dc2883e4 1261
1262 if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ) continue;
1263
d7d7e825 1264 // Check the acceptance for both gammas
01b7fdcc 1265 Bool_t gammaEtaCut = kTRUE;
1266 if(TMath::Abs(daughter0->Eta()) > fV0Reader->GetEtaCut() || TMath::Abs(daughter1->Eta()) > fV0Reader->GetEtaCut() ) gammaEtaCut = kFALSE;
1267 Bool_t gammaRCut = kTRUE;
1268 if(daughter0->R() > fV0Reader->GetMaxRCut() || daughter1->R() > fV0Reader->GetMaxRCut() ) gammaRCut = kFALSE;
d7d7e825 1269
1270
1271
1272 // check for conversions now -> have to pass eta, R and line cut!
1273 Bool_t daughter0Electron = kFALSE;
1274 Bool_t daughter0Positron = kFALSE;
1275 Bool_t daughter1Electron = kFALSE;
1276 Bool_t daughter1Positron = kFALSE;
1277
1278 if(daughter0->GetNDaughters() >= 2){ // first gamma
1279 for(Int_t TrackIndex=daughter0->GetFirstDaughter();TrackIndex<=daughter0->GetLastDaughter();TrackIndex++){
1280 TParticle *tmpDaughter = fStack->Particle(TrackIndex);
1281 if(tmpDaughter->GetUniqueID() == 5){
1282 if(tmpDaughter->GetPdgCode() == 11){
1283 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
1284 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
1285 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
1286 daughter0Electron = kTRUE;
1287 }
1288 }
1289
1290 }
1291 }
1292 else if(tmpDaughter->GetPdgCode() == -11){
1293 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
1294 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
1295 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
1296 daughter0Positron = kTRUE;
1297 }
a68437fb 1298 }
1299 }
d7d7e825 1300 }
1301 }
1302 }
1303 }
1304
1305
1306 if(daughter1->GetNDaughters() >= 2){ // second gamma
1307 for(Int_t TrackIndex=daughter1->GetFirstDaughter();TrackIndex<=daughter1->GetLastDaughter();TrackIndex++){
1308 TParticle *tmpDaughter = fStack->Particle(TrackIndex);
1309 if(tmpDaughter->GetUniqueID() == 5){
1310 if(tmpDaughter->GetPdgCode() == 11){
1311 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
1312 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
1313 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
1314 daughter1Electron = kTRUE;
1315 }
1316 }
1317
1318 }
1319 }
1320 else if(tmpDaughter->GetPdgCode() == -11){
1321 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
1322 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
1323 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
1324 daughter1Positron = kTRUE;
1325 }
1326 }
1327
1328 }
1329
1330 }
1331 }
1332 }
1333 }
1334
a0b94e5c 1335
d7d7e825 1336 if(particle->GetPdgCode()==111){ //Pi0
1337 if( iTracks >= fStack->GetNprimary()){
1338 fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());
1339 fHistograms->FillHistogram("MC_Pi0_Secondaries_Rapid", rapidity);
1340 fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);
1341 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());
1342 fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());
1343 fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());
1344 fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());
1345 fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1346 fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling
1347
01b7fdcc 1348 if(gammaEtaCut && gammaRCut){
d7d7e825 1349 //if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
1350 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1351 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
1352 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1353 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1354 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1355 }
1356 }
1357 }
1358 else{
1359 fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());
1360 fHistograms->FillHistogram("MC_Pi0_Rapid", rapidity);
1361 fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);
1362 fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());
dc2883e4 1363 fHistograms->FillHistogram("MC_Pi0_Pt_vs_Rapid", particle->Pt(),rapidity);
d7d7e825 1364 fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());
1365 fHistograms->FillHistogram("MC_Pi0_R", particle->R());
1366 fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());
1367 fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1368 fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
63e16c52 1369 if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_Fiducial", particle->Pt());
f3026589 1370
1371 if(!fV0Reader->GetIsHeavyIon()) {
1372
1373 switch(GetProcessType(fGCMCEvent)){
1374 case kProcSD:
1375 fHistograms->FillHistogram("MC_SD_EvtQ3_Pi0_Pt", particle->Pt());
1376 break;
1377 case kProcDD:
1378 fHistograms->FillHistogram("MC_DD_EvtQ3_Pi0_Pt", particle->Pt());
1379 break;
1380 case kProcND:
1381 fHistograms->FillHistogram("MC_ND_EvtQ3_Pi0_Pt", particle->Pt());
1382 break;
1383 default:
1384 AliError("Unknown Process");
1385 }
1386 }
01b7fdcc 1387 if(gammaEtaCut && gammaRCut){
d7d7e825 1388 // if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
1389 fHistograms->FillHistogram("MC_Pi0_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1390 fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
63e16c52 1391 if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_withinAcceptance_Fiducial", particle->Pt());
1392
d7d7e825 1393 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1394 fHistograms->FillHistogram("MC_Pi0_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1395 fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1396 fHistograms->FillHistogram("MC_Pi0_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
10e3319b 1397 fHistograms->FillHistogram("MC_Pi0_ConvGamma_OpeningAngle_Pt", particle->Pt(),GetMCOpeningAngle(daughter0,daughter1));
1398 fHistograms->FillHistogram("MC_Pi0_ConvGamma_PtGamma_Pt", particle->Pt(),daughter0->Pt());
1399 fHistograms->FillHistogram("MC_Pi0_ConvGamma_PtGamma_Pt", particle->Pt(),daughter1->Pt());
1400
48682642 1401 Double_t alfa=0.;
77ac6f3e 1402 if((daughter0->Energy()+daughter1->Energy()) > 0.){
48682642 1403 alfa= TMath::Abs((daughter0->Energy()-daughter1->Energy())/(daughter0->Energy()+daughter1->Energy()));
1404 }
1405 fHistograms->FillHistogram("MC_Pi0_alpha",alfa);
1406 if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_ConvGamma_withinAcceptance_Fiducial", particle->Pt());
1407
d7d7e825 1408 }
1409 }
1410 }
1411 }
1412
1413 if(particle->GetPdgCode()==221){ //Eta
1414 fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());
1415 fHistograms->FillHistogram("MC_Eta_Rapid", rapidity);
1416 fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);
1417 fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());
dc2883e4 1418 fHistograms->FillHistogram("MC_Eta_Pt_vs_Rapid", particle->Pt(),rapidity);
d7d7e825 1419 fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());
1420 fHistograms->FillHistogram("MC_Eta_R", particle->R());
1421 fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());
1422 fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1423 fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
1424
01b7fdcc 1425 if(gammaEtaCut && gammaRCut){
a0b94e5c 1426 // if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
d7d7e825 1427 fHistograms->FillHistogram("MC_Eta_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1428 fHistograms->FillHistogram("MC_Eta_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
1429 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1430 fHistograms->FillHistogram("MC_Eta_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1431 fHistograms->FillHistogram("MC_Eta_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1432 fHistograms->FillHistogram("MC_Eta_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
10e3319b 1433 fHistograms->FillHistogram("MC_Eta_ConvGamma_OpeningAngle_Pt", particle->Pt(),GetMCOpeningAngle(daughter0,daughter1));
1434 fHistograms->FillHistogram("MC_Eta_ConvGamma_PtGamma_Pt", particle->Pt(),daughter0->Pt());
1435 fHistograms->FillHistogram("MC_Eta_ConvGamma_PtGamma_Pt", particle->Pt(),daughter1->Pt());
1436
dcdc851f 1437 Double_t alfa=0.;
1438 if((daughter0->Energy()+daughter1->Energy()) > 0.){
1439 alfa= TMath::Abs((daughter0->Energy()-daughter1->Energy())/(daughter0->Energy()+daughter1->Energy()));
1440 }
1441 fHistograms->FillHistogram("MC_Eta_alpha",alfa);
1442
d7d7e825 1443 }
1444
1445 }
1446
1447 }
1448
1449 // all motherparticles with 2 gammas as daughters
1450 fHistograms->FillHistogram("MC_Mother_R", particle->R());
1451 fHistograms->FillHistogram("MC_Mother_ZR", particle->Vz(),particle->R());
1452 fHistograms->FillHistogram("MC_Mother_XY", particle->Vx(),particle->Vy());
1453 fHistograms->FillHistogram("MC_Mother_Mass", particle->GetCalcMass());
1454 fHistograms->FillHistogram("MC_Mother_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1455 fHistograms->FillHistogram("MC_Mother_Energy", particle->Energy());
1456 fHistograms->FillHistogram("MC_Mother_Pt", particle->Pt());
1457 fHistograms->FillHistogram("MC_Mother_Eta", particle->Eta());
1458 fHistograms->FillHistogram("MC_Mother_Rapid", rapidity);
1459 fHistograms->FillHistogram("MC_Mother_Phi",tmpPhi);
1460 fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt",particle->GetMass(),particle->Pt());
a0b94e5c 1461
01b7fdcc 1462 if(gammaEtaCut && gammaRCut){
a0b94e5c 1463 // if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
d7d7e825 1464 fHistograms->FillHistogram("MC_Mother_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1465 fHistograms->FillHistogram("MC_Mother_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
1466 fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_withinAcceptance",particle->GetMass(),particle->Pt());
1467 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1468 fHistograms->FillHistogram("MC_Mother_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1469 fHistograms->FillHistogram("MC_Mother_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1470 fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_ConvGamma_withinAcceptance",particle->GetMass(),particle->Pt());
a0b94e5c 1471
d7d7e825 1472 }
1473
1474
1475 } // end passed R and eta cut
a0b94e5c 1476
d7d7e825 1477 } // end if(particle->GetNDaughters() == 2)
1478
1479 }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)
a0b94e5c 1480
d7d7e825 1481} // end ProcessMCData
1482
1483
1484
1485void AliAnalysisTaskGammaConversion::FillNtuple(){
1486 //Fills the ntuple with the different values
a0b94e5c 1487
d7d7e825 1488 if(fGammaNtuple == NULL){
1489 return;
1490 }
1491 Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
1492 for(Int_t i=0;i<numberOfV0s;i++){
1493 Float_t values[27] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1494 AliESDv0 * cV0 = fV0Reader->GetV0(i);
1495 Double_t negPID=0;
1496 Double_t posPID=0;
1497 fV0Reader->GetPIDProbability(negPID,posPID);
1498 values[0]=cV0->GetOnFlyStatus();
1499 values[1]=fV0Reader->CheckForPrimaryVertex();
1500 values[2]=negPID;
1501 values[3]=posPID;
1502 values[4]=fV0Reader->GetX();
1503 values[5]=fV0Reader->GetY();
1504 values[6]=fV0Reader->GetZ();
1505 values[7]=fV0Reader->GetXYRadius();
1506 values[8]=fV0Reader->GetMotherCandidateNDF();
1507 values[9]=fV0Reader->GetMotherCandidateChi2();
1508 values[10]=fV0Reader->GetMotherCandidateEnergy();
1509 values[11]=fV0Reader->GetMotherCandidateEta();
1510 values[12]=fV0Reader->GetMotherCandidatePt();
1511 values[13]=fV0Reader->GetMotherCandidateMass();
1512 values[14]=fV0Reader->GetMotherCandidateWidth();
1513 // values[15]=fV0Reader->GetMotherMCParticle()->Pt(); MOVED TO THE END, HAS TO BE CALLED AFTER HasSameMother NB: still has the same entry in the array
1514 values[16]=fV0Reader->GetOpeningAngle();
1515 values[17]=fV0Reader->GetNegativeTrackEnergy();
1516 values[18]=fV0Reader->GetNegativeTrackPt();
1517 values[19]=fV0Reader->GetNegativeTrackEta();
1518 values[20]=fV0Reader->GetNegativeTrackPhi();
1519 values[21]=fV0Reader->GetPositiveTrackEnergy();
1520 values[22]=fV0Reader->GetPositiveTrackPt();
1521 values[23]=fV0Reader->GetPositiveTrackEta();
1522 values[24]=fV0Reader->GetPositiveTrackPhi();
1523 values[25]=fV0Reader->HasSameMCMother();
1524 if(values[25] != 0){
1525 values[26]=fV0Reader->GetMotherMCParticlePDGCode();
1526 values[15]=fV0Reader->GetMotherMCParticle()->Pt();
1527 }
1528 fTotalNumberOfAddedNtupleEntries++;
1529 fGammaNtuple->Fill(values);
1530 }
1531 fV0Reader->ResetV0IndexNumber();
1532
1533}
1534
1535void AliAnalysisTaskGammaConversion::ProcessV0sNoCut(){
1536 // Process all the V0's without applying any cuts to it
a0b94e5c 1537
d7d7e825 1538 Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
1539 for(Int_t i=0;i<numberOfV0s;i++){
1540 /*AliESDv0 * cV0 = */fV0Reader->GetV0(i);
a0b94e5c 1541
d7d7e825 1542 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
cb90a330 1543 continue;
d7d7e825 1544 }
9640a3d1 1545
77880bd8 1546 // if( !fV0Reader->GetV0(i)->GetOnFlyStatus()){
1547 if( !fV0Reader->CheckV0FinderStatus(i)){
cb90a330 1548 continue;
9640a3d1 1549 }
1550
cb90a330 1551
1552 if( !((fV0Reader->GetNegativeESDTrack())->GetStatus() & AliESDtrack::kTPCrefit) ||
1553 !((fV0Reader->GetPositiveESDTrack())->GetStatus() & AliESDtrack::kTPCrefit) ){
1554 continue;
9640a3d1 1555 }
cb90a330 1556
ebcfaa7e 1557 if( fV0Reader->GetNegativeESDTrack()->GetSign()== fV0Reader->GetPositiveESDTrack()->GetSign()){
1558 continue;
1559 }
1560
1561 if( fV0Reader->GetNegativeESDTrack()->GetKinkIndex(0) > 0 ||
1562 fV0Reader->GetPositiveESDTrack()->GetKinkIndex(0) > 0) {
1563 continue;
1564 }
9c1cb6f7 1565 if(TMath::Abs(fV0Reader->GetMotherCandidateEta())> fV0Reader->GetEtaCut()){
1566 continue;
1567 }
1568 if(TMath::Abs(fV0Reader->GetPositiveTrackEta())> fV0Reader->GetEtaCut()){
1569 continue;
1570 }
1571 if(TMath::Abs(fV0Reader->GetNegativeTrackEta())> fV0Reader->GetEtaCut()){
1572 continue;
1573 }
1574 if((TMath::Abs(fV0Reader->GetZ())*fV0Reader->GetLineCutZRSlope())-fV0Reader->GetLineCutZValue() > fV0Reader->GetXYRadius() ){ // cuts out regions where we do not reconstruct
1575 continue;
1576 }
f3026589 1577
1578 fHistograms->FillHistogram("ESD_NoCutAllV0_Pt", fV0Reader->GetMotherCandidatePt());
1579
d7d7e825 1580 if(fDoMCTruth){
a0b94e5c 1581
d7d7e825 1582 if(fV0Reader->HasSameMCMother() == kFALSE){
1583 continue;
1584 }
a0b94e5c 1585
d7d7e825 1586 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1587 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
a0b94e5c 1588
d7d7e825 1589 if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1590 continue;
1591 }
26923b22 1592 if(negativeMC->GetPdgCode() == positiveMC->GetPdgCode()){
d7d7e825 1593 continue;
1594 }
a68437fb 1595
26923b22 1596 if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){ // id 5 is conversion
a68437fb 1597 continue;
1598 }
a0b94e5c 1599
d7d7e825 1600 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
a0b94e5c 1601
d7d7e825 1602 fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1603 fHistograms->FillHistogram("ESD_NoCutConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1604 fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1605 fHistograms->FillHistogram("ESD_NoCutConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1606 fHistograms->FillHistogram("ESD_NoCutConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1607 fHistograms->FillHistogram("ESD_NoCutConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1608 fHistograms->FillHistogram("ESD_NoCutConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1609 fHistograms->FillHistogram("ESD_NoCutConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1610 fHistograms->FillHistogram("ESD_NoCutConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1611 fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
a0b94e5c 1612
d7d7e825 1613 fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1614 fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
1615
1616 fHistograms->FillHistogram("ESD_NoCutConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1617 fHistograms->FillHistogram("ESD_NoCutConversion_R", fV0Reader->GetXYRadius());
1618 fHistograms->FillHistogram("ESD_NoCutConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1619 fHistograms->FillHistogram("ESD_NoCutConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
9640a3d1 1620 fHistograms->FillHistogram("ESD_NoCutConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1621 fHistograms->FillHistogram("ESD_NoCutConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1622 fHistograms->FillHistogram("ESD_NoCutConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1623 fHistograms->FillHistogram("ESD_NoCutConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1624
1625 fHistograms->FillHistogram("ESD_NoCutConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1626 fHistograms->FillHistogram("ESD_NoCutConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1627 fHistograms->FillHistogram("ESD_NoCutConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1628 fHistograms->FillHistogram("ESD_NoCutConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1629
d7d7e825 1630 //store MCTruth properties
1631 fHistograms->FillHistogram("ESD_NoCutConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
1632 fHistograms->FillHistogram("ESD_NoCutConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
1633 fHistograms->FillHistogram("ESD_NoCutConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
d7d7e825 1634 }
1635 }
1636 }
1637 fV0Reader->ResetV0IndexNumber();
1638}
1639
1640void AliAnalysisTaskGammaConversion::ProcessV0s(){
1641 // see header file for documentation
037dc2db 1642
1643
d7d7e825 1644 if(fWriteNtuple == kTRUE){
1645 FillNtuple();
1646 }
1647
1648 Int_t nSurvivingV0s=0;
5ce758b0 1649 fV0Reader->ResetNGoodV0s();
d7d7e825 1650 while(fV0Reader->NextV0()){
1651 nSurvivingV0s++;
1652
e60f3265 1653
1654 TVector3 vtxConv(fV0Reader->GetX(),fV0Reader->GetY(), fV0Reader->GetZ());
1655
d7d7e825 1656 //-------------------------- filling v0 information -------------------------------------
1657 fHistograms->FillHistogram("ESD_Conversion_R", fV0Reader->GetXYRadius());
1658 fHistograms->FillHistogram("ESD_Conversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1659 fHistograms->FillHistogram("ESD_Conversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1660 fHistograms->FillHistogram("ESD_Conversion_OpeningAngle", fV0Reader->GetOpeningAngle());
e60f3265 1661
1662 // Specific histograms for beam pipe studies
1663 if( TMath::Abs(fV0Reader->GetZ()) < fV0Reader->GetLineCutZValue() ){
1664 fHistograms->FillHistogram("ESD_Conversion_XY_BeamPipe", fV0Reader->GetX(),fV0Reader->GetY());
1665 fHistograms->FillHistogram("ESD_Conversion_RPhi_BeamPipe", vtxConv.Phi(),fV0Reader->GetXYRadius());
1666 }
1667
1668
d7d7e825 1669 fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());
1670 fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());
1671 fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());
1672 fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());
ebcfaa7e 1673 fHistograms->FillHistogram("ESD_E_nTPCClusters", fV0Reader->GetNegativeTracknTPCClusters());
1674 fHistograms->FillHistogram("ESD_E_nITSClusters", fV0Reader->GetNegativeTracknITSClusters());
40051a3e 1675 Double_t eClsToF= 0;
1676 if(!fV0Reader->GetUseCorrectedTPCClsInfo()){
1677 if(fV0Reader->GetNegativeTracknTPCFClusters()!=0 ){
1678 eClsToF=(Double_t)fV0Reader->GetNegativeTracknTPCClusters()/(Double_t)fV0Reader->GetNegativeTracknTPCFClusters();
1679 }
1680 }else{
1681 eClsToF= fV0Reader->GetNegativeESDTrack()->GetTPCClusterInfo(2,0,fV0Reader->GetFirstTPCRow(fV0Reader->GetXYRadius()));
1682 }
1683 fHistograms->FillHistogram("ESD_E_nTPCClustersToFP", fV0Reader->GetNegativeTrackP(),eClsToF );
1684 fHistograms->FillHistogram("ESD_E_nTPCClustersToFR", fV0Reader->GetXYRadius(),eClsToF );
1685
1686 if(fV0Reader->GetNegativeTracknTPCClusters()!=0 ){
5ce758b0 1687 fHistograms->FillHistogram("ESD_E_TPCchi2", fV0Reader->GetNegativeTrackTPCchi2()/(Double_t)fV0Reader->GetNegativeTracknTPCClusters());
1688 }
1689
1690
1691
d7d7e825 1692 fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());
1693 fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());
1694 fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());
1695 fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());
ebcfaa7e 1696 fHistograms->FillHistogram("ESD_P_nTPCClusters", fV0Reader->GetPositiveTracknTPCClusters());
1697 fHistograms->FillHistogram("ESD_P_nITSClusters", fV0Reader->GetPositiveTracknITSClusters());
40051a3e 1698 Double_t pClsToF= 0;
1699 if(!fV0Reader->GetUseCorrectedTPCClsInfo()){
1700 if(fV0Reader->GetPositiveTracknTPCFClusters()!=0){
1701 pClsToF = (Double_t)fV0Reader->GetPositiveTracknTPCClusters()/(Double_t)fV0Reader->GetPositiveTracknTPCFClusters();
1702 }
1703 }else{
1704 pClsToF= fV0Reader->GetPositiveESDTrack()->GetTPCClusterInfo(2,0,fV0Reader->GetFirstTPCRow(fV0Reader->GetXYRadius()));
1705 }
1706
1707 fHistograms->FillHistogram("ESD_P_nTPCClustersToFP",fV0Reader->GetPositiveTrackP(), pClsToF);
1708 fHistograms->FillHistogram("ESD_P_nTPCClustersToFR",fV0Reader->GetXYRadius(), pClsToF);
1709
1710 if(fV0Reader->GetPositiveTracknTPCClusters()!=0){
5ce758b0 1711 fHistograms->FillHistogram("ESD_P_TPCchi2", fV0Reader->GetPositiveTrackTPCchi2()/(Double_t)fV0Reader->GetPositiveTracknTPCClusters());
1712 }
1713
1714
40051a3e 1715
d7d7e825 1716 fHistograms->FillHistogram("ESD_ConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1717 fHistograms->FillHistogram("ESD_ConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1718 fHistograms->FillHistogram("ESD_ConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1719 fHistograms->FillHistogram("ESD_ConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1720 fHistograms->FillHistogram("ESD_ConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1721 fHistograms->FillHistogram("ESD_ConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1722 fHistograms->FillHistogram("ESD_ConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1723 fHistograms->FillHistogram("ESD_ConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1724 fHistograms->FillHistogram("ESD_ConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1725 fHistograms->FillHistogram("ESD_ConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1726
1727 fHistograms->FillHistogram("ESD_ConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1728 fHistograms->FillHistogram("ESD_ConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
9640a3d1 1729
1730 fHistograms->FillHistogram("ESD_ConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1731 fHistograms->FillHistogram("ESD_ConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1732 fHistograms->FillHistogram("ESD_ConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1733 fHistograms->FillHistogram("ESD_ConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1734
1735 fHistograms->FillHistogram("ESD_ConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1736 fHistograms->FillHistogram("ESD_ConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1737 fHistograms->FillHistogram("ESD_ConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1738 fHistograms->FillHistogram("ESD_ConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1739
9c1cb6f7 1740 Double_t negPID=0;
1741 Double_t posPID=0;
1742 fV0Reader->GetPIDProbability(negPID,posPID);
1743 fHistograms->FillHistogram("ESD_ConvGamma_E_EProbP",fV0Reader->GetNegativeTrackP(),negPID);
1744 fHistograms->FillHistogram("ESD_ConvGamma_P_EProbP",fV0Reader->GetPositiveTrackP(),posPID);
1745
1746 Double_t negPIDmupi=0;
1747 Double_t posPIDmupi=0;
1748 fV0Reader->GetPIDProbabilityMuonPion(negPIDmupi,posPIDmupi);
1749 fHistograms->FillHistogram("ESD_ConvGamma_E_mupiProbP",fV0Reader->GetNegativeTrackP(),negPIDmupi);
1750 fHistograms->FillHistogram("ESD_ConvGamma_P_mupiProbP",fV0Reader->GetPositiveTrackP(),posPIDmupi);
1751
48682642 1752 Double_t armenterosQtAlfa[2];
1753 fV0Reader->GetArmenterosQtAlfa(fV0Reader-> GetNegativeKFParticle(),
1754 fV0Reader-> GetPositiveKFParticle(),
1755 fV0Reader->GetMotherCandidateKFCombination(),
1756 armenterosQtAlfa);
1757
1758 fHistograms->FillHistogram("ESD_ConvGamma_alfa_qt",armenterosQtAlfa[1],armenterosQtAlfa[0]);
1759
037dc2db 1760
d7d7e825 1761 // begin mapping
1762 Int_t rBin = fHistograms->GetRBin(fV0Reader->GetXYRadius());
9640a3d1 1763 Int_t zBin = fHistograms->GetZBin(fV0Reader->GetZ());
d7d7e825 1764 Int_t phiBin = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());
9c1cb6f7 1765 Double_t rFMD=30;
64f4118c 1766 Double_t rITSTPCMin=50;
1767 Double_t rITSTPCMax=80;
e60f3265 1768
9640a3d1 1769
ebcfaa7e 1770 // Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();
d7d7e825 1771
1772 TString nameESDMappingPhiR="";
e158cbc3 1773 nameESDMappingPhiR.Form("ESD_Conversion_Mapping_Phi%02d_R%02d",phiBin,rBin);
ebcfaa7e 1774 //fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);
d7d7e825 1775
1776 TString nameESDMappingPhi="";
e158cbc3 1777 nameESDMappingPhi.Form("ESD_Conversion_Mapping_Phi%02d",phiBin);
ebcfaa7e 1778 //fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);
d7d7e825 1779
1780 TString nameESDMappingR="";
e158cbc3 1781 nameESDMappingR.Form("ESD_Conversion_Mapping_R%02d",rBin);
ebcfaa7e 1782 //fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);
d7d7e825 1783
1784 TString nameESDMappingPhiInR="";
e158cbc3 1785 nameESDMappingPhiInR.Form("ESD_Conversion_Mapping_Phi_in_R_%02d",rBin);
9640a3d1 1786 // fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
1787 fHistograms->FillHistogram(nameESDMappingPhiInR, vtxConv.Phi());
1788
1789 TString nameESDMappingZInR="";
1790 nameESDMappingZInR.Form("ESD_Conversion_Mapping_Z_in_R_%02d",rBin);
1791 fHistograms->FillHistogram(nameESDMappingZInR, fV0Reader->GetZ());
1792
1793 TString nameESDMappingPhiInZ="";
1794 nameESDMappingPhiInZ.Form("ESD_Conversion_Mapping_Phi_in_Z_%02d",zBin);
1795 // fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
1796 fHistograms->FillHistogram(nameESDMappingPhiInZ, vtxConv.Phi());
1797
9c1cb6f7 1798 if(fV0Reader->GetXYRadius()<rFMD){
1799 TString nameESDMappingFMDPhiInZ="";
1800 nameESDMappingFMDPhiInZ.Form("ESD_Conversion_Mapping_FMD_Phi_in_Z_%02d",zBin);
1801 fHistograms->FillHistogram(nameESDMappingFMDPhiInZ, vtxConv.Phi());
1802 }
1803
64f4118c 1804 if(fV0Reader->GetXYRadius()>rITSTPCMin && fV0Reader->GetXYRadius()<rITSTPCMax){
1805 TString nameESDMappingITSTPCPhiInZ="";
1806 nameESDMappingITSTPCPhiInZ.Form("ESD_Conversion_Mapping_ITSTPC_Phi_in_Z_%02d",zBin);
1807 fHistograms->FillHistogram(nameESDMappingITSTPCPhiInZ, vtxConv.Phi());
1808 }
9c1cb6f7 1809
9640a3d1 1810 TString nameESDMappingRInZ="";
1811 nameESDMappingRInZ.Form("ESD_Conversion_Mapping_R_in_Z_%02d",zBin);
1812 fHistograms->FillHistogram(nameESDMappingRInZ, fV0Reader->GetXYRadius());
1813
1814 if(fV0Reader->GetMotherCandidatePt() > fLowPtMapping && fV0Reader->GetMotherCandidatePt()< fHighPtMapping){
1815 TString nameESDMappingMidPtPhiInR="";
1816 nameESDMappingMidPtPhiInR.Form("ESD_Conversion_Mapping_MidPt_Phi_in_R_%02d",rBin);
1817 fHistograms->FillHistogram(nameESDMappingMidPtPhiInR, vtxConv.Phi());
1818
1819 TString nameESDMappingMidPtZInR="";
1820 nameESDMappingMidPtZInR.Form("ESD_Conversion_Mapping_MidPt_Z_in_R_%02d",rBin);
1821 fHistograms->FillHistogram(nameESDMappingMidPtZInR, fV0Reader->GetZ());
1822
1823 TString nameESDMappingMidPtPhiInZ="";
1824 nameESDMappingMidPtPhiInZ.Form("ESD_Conversion_Mapping_MidPt_Phi_in_Z_%02d",zBin);
1825 fHistograms->FillHistogram(nameESDMappingMidPtPhiInZ, vtxConv.Phi());
9c1cb6f7 1826 if(fV0Reader->GetXYRadius()<rFMD){
1827 TString nameESDMappingMidPtFMDPhiInZ="";
1828 nameESDMappingMidPtFMDPhiInZ.Form("ESD_Conversion_Mapping_MidPt_FMD_Phi_in_Z_%02d",zBin);
1829 fHistograms->FillHistogram(nameESDMappingMidPtFMDPhiInZ, vtxConv.Phi());
1830 }
1831
9640a3d1 1832
1833 TString nameESDMappingMidPtRInZ="";
1834 nameESDMappingMidPtRInZ.Form("ESD_Conversion_Mapping_MidPt_R_in_Z_%02d",zBin);
1835 fHistograms->FillHistogram(nameESDMappingMidPtRInZ, fV0Reader->GetXYRadius());
1836
1837 }
1838
1839
d7d7e825 1840 // end mapping
1841
6c84d371 1842 new((*fKFReconstructedGammasTClone)[fKFReconstructedGammasTClone->GetEntriesFast()]) AliKFParticle(*fV0Reader->GetMotherCandidateKFCombination());
037dc2db 1843 fKFReconstructedGammasV0Index.push_back(fV0Reader->GetCurrentV0IndexNumber()-1);
6c84d371 1844 // fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());
d7d7e825 1845 fElectronv1.push_back(fV0Reader->GetCurrentV0()->GetPindex());
1846 fElectronv2.push_back(fV0Reader->GetCurrentV0()->GetNindex());
a0b94e5c 1847
d7d7e825 1848 //----------------------------------- checking for "real" conversions (MC match) --------------------------------------
1849 if(fDoMCTruth){
f3026589 1850 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1851 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1852
1853 if(fV0Reader->HasSameMCMother() == kFALSE){
1854 fHistograms->FillHistogram("ESD_TrueConvCombinatorial_R", fV0Reader->GetXYRadius());
1855 fHistograms->FillHistogram("ESD_TrueConvCombinatorial_Pt", fV0Reader->GetMotherCandidatePt());
1856 fHistograms->FillHistogram("ESD_TrueConvCombinatorialDaughter_Pt", negativeMC->Pt(),positiveMC->Pt());
1857 if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
1858 fHistograms->FillHistogram("ESD_TrueConvCombinatorialElec_R", fV0Reader->GetXYRadius());
1859 fHistograms->FillHistogram("ESD_TrueConvCombinatorialElec_Pt", fV0Reader->GetMotherCandidatePt());
1860 }
1861 if(TMath::Abs(negativeMC->GetPdgCode())==211 && TMath::Abs(positiveMC->GetPdgCode())==211){
1862 fHistograms->FillHistogram("ESD_TrueConvCombinatorialPi_R", fV0Reader->GetXYRadius());
1863 fHistograms->FillHistogram("ESD_TrueConvCombinatorialPi_Pt", fV0Reader->GetMotherCandidatePt());
1864 fHistograms->FillHistogram("ESD_TrueConvCombinatorialPiDaughter_Pt", negativeMC->Pt(),positiveMC->Pt());
1865 }
1866 if((TMath::Abs(negativeMC->GetPdgCode())==211 && TMath::Abs(positiveMC->GetPdgCode())==2211) ||
1867 (TMath::Abs(negativeMC->GetPdgCode())==2212 && TMath::Abs(positiveMC->GetPdgCode())==211)){
1868 fHistograms->FillHistogram("ESD_TrueConvCombinatorialPiP_R", fV0Reader->GetXYRadius());
1869 fHistograms->FillHistogram("ESD_TrueConvCombinatorialPiP_Pt", fV0Reader->GetMotherCandidatePt());
1870 fHistograms->FillHistogram("ESD_TrueConvCombinatorialPiPDaughter_Pt", negativeMC->Pt(),positiveMC->Pt());
1871 }
1872 if((TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==211) ||
1873 (TMath::Abs(negativeMC->GetPdgCode())==211 && TMath::Abs(positiveMC->GetPdgCode())==11)){
1874 fHistograms->FillHistogram("ESD_TrueConvCombinatorialElecPi_R", fV0Reader->GetXYRadius());
1875 fHistograms->FillHistogram("ESD_TrueConvCombinatorialElecPi_Pt", fV0Reader->GetMotherCandidatePt());
1876 }
1877 continue;
1878 }
1879 // Moved up to check true electron background
1880 // TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1881 // TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1882
1883 if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1884 fHistograms->FillHistogram("ESD_TrueConvHadronicBck_R", fV0Reader->GetXYRadius());
1885 fHistograms->FillHistogram("ESD_TrueConvHadronicBck_Pt", fV0Reader->GetMotherCandidatePt());
1886 fHistograms->FillHistogram("ESD_TrueConvHadronicBckDaughter_Pt", negativeMC->Pt(),positiveMC->Pt());
1887 if((TMath::Abs(negativeMC->GetPdgCode())==211 && TMath::Abs(positiveMC->GetPdgCode())==2211) ||
1888 (TMath::Abs(negativeMC->GetPdgCode())==2212 && TMath::Abs(positiveMC->GetPdgCode())==211)){
1889 fHistograms->FillHistogram("ESD_TrueConvLambda_R", fV0Reader->GetXYRadius());
1890 fHistograms->FillHistogram("ESD_TrueConvLambda_Pt", fV0Reader->GetMotherCandidatePt());
1891 }
1892 if(TMath::Abs(negativeMC->GetPdgCode())==211 && TMath::Abs(positiveMC->GetPdgCode())==211){
1893 fHistograms->FillHistogram("ESD_TrueConvMeson_R", fV0Reader->GetXYRadius());
1894 fHistograms->FillHistogram("ESD_TrueConvMeson_Pt", fV0Reader->GetMotherCandidatePt());
1895 }
1896 continue;
1897 }
1898
a68437fb 1899
f3026589 1900 if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
1901 continue;
1902 }
8de8695b 1903
1904 UInt_t statusPos = fV0Reader->GetPositiveESDTrack()->GetStatus();
1905 UInt_t statusNeg = fV0Reader->GetNegativeESDTrack()->GetStatus();
1906 UChar_t itsPixelPos = fV0Reader->GetPositiveESDTrack()->GetITSClusterMap();
1907 UChar_t itsPixelNeg = fV0Reader->GetNegativeESDTrack()->GetITSClusterMap();
1908
1909 // Using the UniqueID Phojet does not get the Dalitz right
1910 // if( (negativeMC->GetUniqueID() == 4 && positiveMC->GetUniqueID() ==4) ||
1911 // (negativeMC->GetUniqueID() == 0 && positiveMC->GetUniqueID() ==0) ){// fill r distribution for Dalitz decays
1912 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 111){ //pi0
1913 fHistograms->FillHistogram("ESD_TrueDalitzContamination_R", fV0Reader->GetXYRadius());
1914 fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_R", fV0Reader->GetXYRadius());
1915 //--------Histos for HFE
1916
1917 if(statusPos & AliESDtrack::kTOFpid){
1918 fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_SinglePos_R", fV0Reader->GetXYRadius());
1919 if( TESTBIT(itsPixelPos, 0) ){
1920 fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_SinglePos_kFirst_R", fV0Reader->GetXYRadius());
1921 }
037dc2db 1922 }
8de8695b 1923 if(statusNeg & AliESDtrack::kTOFpid){
1924 fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_SingleNeg_R", fV0Reader->GetXYRadius());
1925 if( TESTBIT(itsPixelNeg, 0) ){
1926 fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_SingleNeg_kFirst_R", fV0Reader->GetXYRadius());
1927 }
64f4118c 1928 }
8de8695b 1929 //--------------------------------------------------------
64f4118c 1930
037dc2db 1931 }
8de8695b 1932 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 221){ //eta
1933 fHistograms->FillHistogram("ESD_TrueConvDalitzEta_R", fV0Reader->GetXYRadius());
1934 }
1935
1936 //}
a68437fb 1937
1938 if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){// check if the daughters come from a conversion
1939 continue;
1940 }
a0b94e5c 1941
d7d7e825 1942 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
a68437fb 1943
26923b22 1944 if(fDoCF){
1945 Double_t containerInput[3];
1946 containerInput[0] = fV0Reader->GetMotherCandidatePt();
1947 containerInput[1] = fV0Reader->GetMotherCandidateEta();
1948 containerInput[2] = fV0Reader->GetMotherCandidateMass();
1949 fCFManager->GetParticleContainer()->Fill(containerInput,kStepTrueGamma); // for CF
1950 }
1951
d7d7e825 1952 fHistograms->FillHistogram("ESD_TrueConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1953 fHistograms->FillHistogram("ESD_TrueConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1954 fHistograms->FillHistogram("ESD_TrueConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1955 fHistograms->FillHistogram("ESD_TrueConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1956 fHistograms->FillHistogram("ESD_TrueConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1957 fHistograms->FillHistogram("ESD_TrueConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1958 fHistograms->FillHistogram("ESD_TrueConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1959 fHistograms->FillHistogram("ESD_TrueConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1960 fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1961 fHistograms->FillHistogram("ESD_TrueConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1962 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters());
1963 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters());
1964 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters(),fV0Reader->GetMotherCandidateMass());
1965 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters(),fV0Reader->GetMotherCandidateMass());
a0b94e5c 1966
d7d7e825 1967 fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1968 fHistograms->FillHistogram("ESD_TrueConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
a0b94e5c 1969
40051a3e 1970 fHistograms->FillHistogram("ESD_TrueConversion_E_nTPCClustersToFR", fV0Reader->GetXYRadius(),eClsToF );
1971 fHistograms->FillHistogram("ESD_TrueConversion_P_nTPCClustersToFR",fV0Reader->GetXYRadius(), pClsToF);
86e91247 1972
d7d7e825 1973 fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1974 fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());
1975 fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1976 fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
9640a3d1 1977
8de8695b 1978 //----Histos for HFE--------------------------------------
1979
1980 if(statusPos & AliESDtrack::kTOFpid){
1981 fHistograms->FillHistogram("ESD_TrueConversion_SinglePos_R", positiveMC->R(),fV0Reader->GetPositiveMCParticle()->Pt());
1982 if( TESTBIT(itsPixelPos, 0) ){
1983 fHistograms->FillHistogram("ESD_TrueConversion_SinglePos_kFirst_R", positiveMC->R(),fV0Reader->GetPositiveMCParticle()->Pt());
1984 }
1985 }
1986 if(statusNeg & AliESDtrack::kTOFpid){
1987 fHistograms->FillHistogram("ESD_TrueConversion_SingleNeg_R", negativeMC->R(),fV0Reader->GetNegativeMCParticle()->Pt());
1988 if( TESTBIT(itsPixelNeg, 0) ){
1989 fHistograms->FillHistogram("ESD_TrueConversion_SingleNeg_kFirst_R", negativeMC->R(),fV0Reader->GetNegativeMCParticle()->Pt());
1990 }
1991 }
1992 //--------------------------------------------------------
1993
9640a3d1 1994 fHistograms->FillHistogram("ESD_TrueConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1995 fHistograms->FillHistogram("ESD_TrueConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1996 fHistograms->FillHistogram("ESD_TrueConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1997 fHistograms->FillHistogram("ESD_TrueConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
96ade8ef 1998 if (fV0Reader->GetMotherCandidateP() != 0) {
1999 fHistograms->FillHistogram("ESD_TrueConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
2000 fHistograms->FillHistogram("ESD_TrueConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
2001 } else { cout << "Error::fV0Reader->GetNegativeTrackP() == 0 !!!" << endl; }
2002
9640a3d1 2003 fHistograms->FillHistogram("ESD_TrueConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
2004 fHistograms->FillHistogram("ESD_TrueConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
70ef88b5 2005 fHistograms->FillHistogram("ESD_TrueConvGamma_alfa_qt",armenterosQtAlfa[1],armenterosQtAlfa[0]);
2006
9640a3d1 2007
a0b94e5c 2008
d7d7e825 2009 //store MCTruth properties
2010 fHistograms->FillHistogram("ESD_TrueConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
2011 fHistograms->FillHistogram("ESD_TrueConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
2012 fHistograms->FillHistogram("ESD_TrueConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
48682642 2013
d7d7e825 2014 //resolution
2015 Double_t mcpt = fV0Reader->GetMotherMCParticle()->Pt();
2016 Double_t esdpt = fV0Reader->GetMotherCandidatePt();
9c1cb6f7 2017 Double_t resdPt = 0.;
4a6157dc 2018 if(mcpt > 0){
9c1cb6f7 2019 resdPt = ((esdpt - mcpt)/mcpt)*100.;
d7d7e825 2020 }
4a6157dc 2021 else if(mcpt < 0){
2022 cout<<"Pt of MC particle is negative, this will cause wrong calculation of resPt"<<endl;
2023 }
d7d7e825 2024
ca6d4600 2025 fHistograms->FillHistogram("Resolution_Gamma_dPt_Pt", mcpt, resdPt);
d7d7e825 2026 fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);
2027 fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);
ca6d4600 2028 fHistograms->FillHistogram("Resolution_Gamma_dPt_Phi", fV0Reader->GetMotherCandidatePhi(), resdPt);
d7d7e825 2029
9c1cb6f7 2030 Double_t resdZ = 0.;
d7d7e825 2031 if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){
9c1cb6f7 2032 resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100.;
d7d7e825 2033 }
9c1cb6f7 2034 Double_t resdZAbs = 0.;
037dc2db 2035 resdZAbs = (fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz());
2036
2037 fHistograms->FillHistogram("Resolution_dZAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdZAbs);
d7d7e825 2038 fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);
2039 fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());
2040 fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());
ca6d4600 2041
2042 // new for dPt_Pt-histograms for Electron and Positron
2043 Double_t mcEpt = fV0Reader->GetNegativeMCParticle()->Pt();
9c1cb6f7 2044 Double_t resEdPt = 0.;
ca6d4600 2045 if (mcEpt > 0){
9c1cb6f7 2046 resEdPt = ((fV0Reader->GetNegativeTrackPt()-mcEpt)/mcEpt)*100.;
ca6d4600 2047 }
2048 UInt_t statusN = fV0Reader->GetNegativeESDTrack()->GetStatus();
2049 // AliESDtrack * negTrk = fV0Reader->GetNegativeESDTrack();
2050 UInt_t kTRDoutN = (statusN & AliESDtrack::kTRDout);
2051
77ac6f3e 2052 Int_t nITSclsE= fV0Reader->GetNegativeTracknITSClusters();
ca6d4600 2053 // filling Resolution_Pt_dPt with respect to the Number of ITS clusters for Positrons
77ac6f3e 2054 switch(nITSclsE){
ca6d4600 2055 case 0: // 0 ITS clusters
2056 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS0", mcEpt, resEdPt);
2057 break;
2058 case 1: // 1 ITS cluster
2059 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS1", mcEpt, resEdPt);
2060 break;
2061 case 2: // 2 ITS clusters
2062 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS2", mcEpt, resEdPt);
2063 break;
2064 case 3: // 3 ITS clusters
2065 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS3", mcEpt, resEdPt);
2066 break;
2067 case 4: // 4 ITS clusters
2068 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS4", mcEpt, resEdPt);
2069 break;
2070 case 5: // 5 ITS clusters
2071 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS5", mcEpt, resEdPt);
2072 break;
2073 case 6: // 6 ITS clusters
2074 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS6", mcEpt, resEdPt);
2075 break;
2076 }
2077 //Filling histograms with respect to Electron resolution
2078 fHistograms->FillHistogram("Resolution_E_dPt_Pt", mcEpt, resEdPt);
2079 fHistograms->FillHistogram("Resolution_E_dPt_Phi", fV0Reader->GetNegativeTrackPhi(), resEdPt);
2080 if(kTRDoutN){
2081 fHistograms->FillHistogram("Resolution_E_nTRDtracklets_ESDPt", fV0Reader->GetNegativeTrackPt(), fV0Reader->GetNegativeESDTrack()->GetTRDntracklets());
2082 fHistograms->FillHistogram("Resolution_E_nTRDtracklets_MCPt", mcEpt, fV0Reader->GetNegativeESDTrack()->GetTRDntracklets());
2083 fHistograms->FillHistogram("Resolution_E_nTRDclusters_ESDPt",fV0Reader->GetNegativeTrackPt(), fV0Reader->GetNegativeESDTrack()->GetTRDncls());
2084 fHistograms->FillHistogram("Resolution_E_nTRDclusters_MCPt",mcEpt, fV0Reader->GetNegativeESDTrack()->GetTRDncls());
2085 fHistograms->FillHistogram("Resolution_E_TRDsignal_ESDPt", fV0Reader->GetNegativeTrackPt(), fV0Reader->GetNegativeESDTrack()->GetTRDsignal());
2086 }
2087
2088 Double_t mcPpt = fV0Reader->GetPositiveMCParticle()->Pt();
2089 Double_t resPdPt = 0;
2090 if (mcPpt > 0){
9c1cb6f7 2091 resPdPt = ((fV0Reader->GetPositiveTrackPt()-mcPpt)/mcPpt)*100.;
ca6d4600 2092 }
2093
2094 UInt_t statusP = fV0Reader->GetPositiveESDTrack()->GetStatus();
2095// AliESDtrack * posTr= fV0Reader->GetPositiveESDTrack();
2096 UInt_t kTRDoutP = (statusP & AliESDtrack::kTRDout);
2097
77ac6f3e 2098 Int_t nITSclsP = fV0Reader->GetPositiveTracknITSClusters();
ca6d4600 2099 // filling Resolution_Pt_dPt with respect to the Number of ITS clusters for Positrons
77ac6f3e 2100 switch(nITSclsP){
ca6d4600 2101 case 0: // 0 ITS clusters
2102 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS0", mcPpt, resPdPt);
2103 break;
2104 case 1: // 1 ITS cluster
2105 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS1", mcPpt, resPdPt);
2106 break;
2107 case 2: // 2 ITS clusters
2108 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS2", mcPpt, resPdPt);
2109 break;
2110 case 3: // 3 ITS clusters
2111 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS3", mcPpt, resPdPt);
2112 break;
2113 case 4: // 4 ITS clusters
2114 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS4", mcPpt, resPdPt);
2115 break;
2116 case 5: // 5 ITS clusters
2117 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS5", mcPpt, resPdPt);
2118 break;
2119 case 6: // 6 ITS clusters
2120 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS6", mcPpt, resPdPt);
2121 break;
2122 }
2123 //Filling histograms with respect to Positron resolution
2124 fHistograms->FillHistogram("Resolution_P_dPt_Pt", mcPpt, resPdPt);
2125 fHistograms->FillHistogram("Resolution_P_dPt_Phi", fV0Reader->GetPositiveTrackPhi(), resPdPt);
2126 if(kTRDoutP){
2127 fHistograms->FillHistogram("Resolution_P_nTRDtracklets_ESDPt", fV0Reader->GetPositiveTrackPt(), fV0Reader->GetPositiveESDTrack()->GetTRDntracklets());
2128 fHistograms->FillHistogram("Resolution_P_nTRDtracklets_MCPt", mcPpt, fV0Reader->GetPositiveESDTrack()->GetTRDntracklets());
2129 fHistograms->FillHistogram("Resolution_P_nTRDclusters_ESDPt",fV0Reader->GetPositiveTrackPt(), fV0Reader->GetPositiveESDTrack()->GetTRDncls());
2130 fHistograms->FillHistogram("Resolution_P_nTRDclusters_MCPt",mcPpt, fV0Reader->GetPositiveESDTrack()->GetTRDncls());
2131 fHistograms->FillHistogram("Resolution_P_TRDsignal_ESDPt", fV0Reader->GetPositiveTrackPt(), fV0Reader->GetPositiveESDTrack()->GetTRDsignal());
2132 }
2133
2134
9c1cb6f7 2135 Double_t resdR = 0.;
d7d7e825 2136 if(fV0Reader->GetNegativeMCParticle()->R() != 0){
9c1cb6f7 2137 resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100.;
d7d7e825 2138 }
9c1cb6f7 2139 Double_t resdRAbs = 0.;
037dc2db 2140 resdRAbs = (fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R());
2141
2142 fHistograms->FillHistogram("Resolution_dRAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdRAbs);
d7d7e825 2143 fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);
2144 fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());
2145 fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());
ca6d4600 2146 fHistograms->FillHistogram("Resolution_R_dPt", fV0Reader->GetNegativeMCParticle()->R(), resdPt);
2147
9c1cb6f7 2148 Double_t resdPhiAbs=0.;
2149 resdPhiAbs=0.;
037dc2db 2150 resdPhiAbs= (fV0Reader->GetMotherCandidatePhi()-fV0Reader->GetNegativeMCParticle()->Phi());
2151 fHistograms->FillHistogram("Resolution_dPhiAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdPhiAbs);
2152
d7d7e825 2153 }//if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22)
d7d7e825 2154 }//if(fDoMCTruth)
2155 }//while(fV0Reader->NextV0)
2156 fHistograms->FillHistogram("ESD_NumberOfSurvivingV0s", nSurvivingV0s);
2157 fHistograms->FillHistogram("ESD_NumberOfV0s", fV0Reader->GetNumberOfV0s());
b5832f95 2158 fHistograms->FillHistogram("ESD_NumberOfContributorsVtx", fV0Reader->GetNumberOfContributorsVtx());
cb90a330 2159
2160 fV0Reader->ResetV0IndexNumber();
d7d7e825 2161}
2162
6441e967 2163void AliAnalysisTaskGammaConversion::AddToAODBranch(TClonesArray * branch, AliAODPWG4Particle & particle) {
2164 //See header file for documentation
2165
2166 Int_t i = branch->GetEntriesFast();
2167 if(! (fOutputAODClassName.Contains("Correlation")) ) {
2168 new((*branch)[i]) AliAODPWG4Particle(particle);
2169 } else {
2170 new((*branch)[i]) AliAODPWG4ParticleCorrelation(particle);
2171 }
2172}
2173
2174void AliAnalysisTaskGammaConversion::AddToAODBranch(TClonesArray * branch, AliGammaConversionAODObject & particle) {
2175 //See header file for documentation
2176
2177 Int_t i = branch->GetEntriesFast();
2178 new((*branch)[i]) AliGammaConversionAODObject(particle);
2179}
2180
2181void AliAnalysisTaskGammaConversion::AddToAODBranch(TClonesArray * branch, AliAODConversionParticle & particle) {
2182 //See header file for documentation
2183
2184 Int_t i = branch->GetEntriesFast();
2185 new((*branch)[i]) AliAODConversionParticle(particle);
2186}
2187
2188
d7d7e825 2189void AliAnalysisTaskGammaConversion::FillAODWithConversionGammas(){
2190 // Fill AOD with reconstructed Gamma
a0b94e5c 2191
6c84d371 2192 for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
332f1f44 2193 AliKFParticle * gammakf = dynamic_cast<AliKFParticle*>(fKFReconstructedGammasTClone->At(gammaIndex));
6441e967 2194
2195 if(fOutputAODClassName.Contains("AliAODPWG4Particle")) {
2196 AliAODPWG4Particle gamma = AliAODPWG4Particle(gammakf->Px(),gammakf->Py(),gammakf->Pz(), gammakf->E());
2197 //gamma.SetLabel(-1);//How to get the MC label of the reconstructed gamma?
2198 gamma.SetTrackLabel( fElectronv1[gammaIndex], fElectronv2[gammaIndex] ); //How to get the MC label of the 2 electrons that form the gamma?
2199 gamma.SetDetector("CTS"); //tag the gamma as reconstructed in the central barrel
2200 gamma.SetPdg(AliPID::kEleCon); //photon id
2201 gamma.SetTag(-1); //Here I usually put a flag saying that montecarlo says it is prompt, decay fragmentation photon, or hadrons or whatever
5510ebe1 2202 //PH gamma.SetChi2(gammakf->Chi2());
6441e967 2203
2204 AddToAODBranch(fAODGamma, gamma);
2205
2206 } else if(fOutputAODClassName.Contains("ConversionParticle")) {
2207 TLorentzVector momentum(gammakf->Px(),gammakf->Py(),gammakf->Pz(), gammakf->E());
2208 AliAODConversionParticle gamma = AliAODConversionParticle(momentum);
2209 //gamma.SetLabel(-1);//How to get the MC label of the reconstructed gamma?
2210 gamma.SetTrackLabels( fElectronv1[gammaIndex], fElectronv2[gammaIndex] ); //How to get the MC label of the 2 electrons that form the gamma?
2211 //gamma.SetPdg(AliPID::kEleCon); //photon id
2212 //gamma.SetTag(-1); //Here I usually put a flag saying that montecarlo says it is prompt, decay fragmentation photon, or hadrons or whatever
c8b45866 2213 gamma.SetChi2(gammakf->Chi2());
6441e967 2214 gamma.SetTrackLabels( fElectronv1[gammaIndex], fElectronv2[gammaIndex] );
2215 gamma.SetESDEvent(dynamic_cast<AliESDEvent*>(InputEvent()));
2216 AddToAODBranch(fAODGamma, gamma);
2217
2218
2219
2220 } else {
2221 AliGammaConversionAODObject gamma;
2222 gamma.SetPx(gammakf->GetPx());
2223 gamma.SetPy(gammakf->GetPy());
2224 gamma.SetPz(gammakf->GetPz());
e0b4c21c 2225 gamma.SetE(gammakf->GetE());
6441e967 2226 gamma.SetLabel1(fElectronv1[gammaIndex]);
2227 gamma.SetLabel2(fElectronv2[gammaIndex]);
2228 gamma.SetChi2(gammakf->Chi2());
2229 gamma.SetE(gammakf->E());
2230 gamma.SetESDEvent(dynamic_cast<AliESDEvent*>(InputEvent()));
2231 AddToAODBranch(fAODGamma, gamma);
2232 }
a0b94e5c 2233
6441e967 2234 }
d7d7e825 2235}
6272370b 2236void AliAnalysisTaskGammaConversion::ProcessGammasForOmegaMesonAnalysis(){
2237 // omega meson analysis pi0+gamma decay
2238 for(Int_t firstPi0Index=0;firstPi0Index<fKFReconstructedPi0sTClone->GetEntriesFast();firstPi0Index++){
48682642 2239 AliKFParticle * omegaCandidatePi0Daughter = (AliKFParticle *)fKFReconstructedPi0sTClone->At(firstPi0Index);
6272370b 2240 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
48682642 2241
6272370b 2242 AliKFParticle * omegaCandidateGammaDaughter = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2243 if(fGammav1[firstPi0Index]==firstGammaIndex || fGammav2[firstPi0Index]==firstGammaIndex){
2244 continue;
2245 }
2246
48682642 2247 AliKFParticle omegaCandidate(*omegaCandidatePi0Daughter,*omegaCandidateGammaDaughter);
6272370b 2248 Double_t massOmegaCandidate = 0.;
2249 Double_t widthOmegaCandidate = 0.;
2250
48682642 2251 omegaCandidate.GetMass(massOmegaCandidate,widthOmegaCandidate);
6272370b 2252
e9aea39f 2253 if ( massOmegaCandidate > 733 && massOmegaCandidate < 833 ) {
6441e967 2254 //AddOmegaToAOD(&omegaCandidate, massOmegaCandidate, firstPi0Index, firstGammaIndex);
c59360eb 2255 }
2256
48682642 2257 fHistograms->FillHistogram("ESD_Omega_InvMass_vs_Pt",massOmegaCandidate ,omegaCandidate.GetPt());
6272370b 2258 fHistograms->FillHistogram("ESD_Omega_InvMass",massOmegaCandidate);
48682642 2259
2260 //delete omegaCandidate;
2261
2262 }// end of omega reconstruction in pi0+gamma channel
2263
2264 if(fDoJet == kTRUE){
2265 AliKFParticle* negPiKF=NULL;
2266 AliKFParticle* posPiKF=NULL;
2267
2268 // look at the pi+pi+pi0 channel
2269 for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
2270 AliESDtrack* posTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
2271 if (posTrack->GetSign()<0) continue;
c0d9051e 2272 if(TMath::Abs(fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kPion))>2.) continue;
48682642 2273 if (posPiKF) delete posPiKF; posPiKF=NULL;
2274 posPiKF = new AliKFParticle( *(posTrack) ,211);
2275
2276 for(Int_t jCh=0;jCh<fChargedParticles->GetEntriesFast();jCh++){
2277 AliESDtrack* negTrack = (AliESDtrack*)(fChargedParticles->At(jCh));
2278 if( negTrack->GetSign()>0) continue;
c0d9051e 2279 if(TMath::Abs(fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kPion))>2.) continue;
48682642 2280 if (negPiKF) delete negPiKF; negPiKF=NULL;
2281 negPiKF = new AliKFParticle( *(negTrack) ,-211);
2282 AliKFParticle omegaCandidatePipPinPi0(*omegaCandidatePi0Daughter,*posPiKF,*negPiKF);
2283 Double_t massOmegaCandidatePipPinPi0 = 0.;
2284 Double_t widthOmegaCandidatePipPinPi0 = 0.;
2285
2286 omegaCandidatePipPinPi0.GetMass(massOmegaCandidatePipPinPi0,widthOmegaCandidatePipPinPi0);
ef2e2748 2287
2288 if ( massOmegaCandidatePipPinPi0 > 733 && massOmegaCandidatePipPinPi0 < 833 ) {
6441e967 2289 // AddOmegaToAOD(&omegaCandidatePipPinPi0, massOmegaCandidatePipPinPi0, -1, -1);
ef2e2748 2290 }
48682642 2291
2292 fHistograms->FillHistogram("ESD_OmegaPipPinPi0_InvMass_vs_Pt",massOmegaCandidatePipPinPi0 ,omegaCandidatePipPinPi0.GetPt());
2293 fHistograms->FillHistogram("ESD_OmegaPipPinPi0_InvMass",massOmegaCandidatePipPinPi0);
2294
2295 // delete omegaCandidatePipPinPi0;
2296 }
2297 }
48682642 2298
45936300 2299 if (posPiKF) delete posPiKF; posPiKF=NULL; if (negPiKF) delete negPiKF; negPiKF=NULL;
48682642 2300
45936300 2301 } // checking ig gammajet because in that case the chargedparticle list is created
6272370b 2302
6272370b 2303 }
48682642 2304
2305 if(fCalculateBackground){
5ce758b0 2306
2307 AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
2308
2309 Int_t zbin= bgHandler->GetZBinIndex(fV0Reader->GetVertexZ());
2310 Int_t mbin = 0;
2311 if(fUseTrackMultiplicityForBG == kTRUE){
2312 mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2313 }
2314 else{
2315 mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->GetNGoodV0s());
2316 }
2317
2318 AliGammaConversionBGHandler::GammaConversionVertex *bgEventVertex = NULL;
2319
48682642 2320 // Background calculation for the omega
2321 for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
5ce758b0 2322 AliGammaConversionKFVector * previousEventV0s = bgHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);
2323
2324 if(fMoveParticleAccordingToVertex == kTRUE){
2325 bgEventVertex = bgHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
2326 }
48682642 2327 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2328 AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
5ce758b0 2329
2330 if(fMoveParticleAccordingToVertex == kTRUE){
2331 MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
2332 }
2333
48682642 2334 for(Int_t firstPi0Index=0;firstPi0Index<fKFReconstructedPi0sTClone->GetEntriesFast();firstPi0Index++){
2335 AliKFParticle * omegaCandidatePi0Daughter = (AliKFParticle *)fKFReconstructedPi0sTClone->At(firstPi0Index);
2336 AliKFParticle * omegaBckCandidate = new AliKFParticle(*omegaCandidatePi0Daughter,previousGoodV0);
2337 Double_t massOmegaBckCandidate = 0.;
2338 Double_t widthOmegaBckCandidate = 0.;
2339
2340 omegaBckCandidate->GetMass(massOmegaBckCandidate,widthOmegaBckCandidate);
ef2e2748 2341
2342
48682642 2343 fHistograms->FillHistogram("ESD_Omega_Bck_InvMass_vs_Pt",massOmegaBckCandidate ,omegaBckCandidate->GetPt());
2344 fHistograms->FillHistogram("ESD_Omega_Bck_InvMass",massOmegaBckCandidate);
2345
2346 delete omegaBckCandidate;
2347 }
2348 }
2349 }
2350 } // end of checking if background calculation is available
6272370b 2351}
d7d7e825 2352
04bf4381 2353
77ac6f3e 2354void AliAnalysisTaskGammaConversion::AddOmegaToAOD(const AliKFParticle * const omegakf, Double_t mass, Int_t omegaDaughter, Int_t gammaDaughter) {
04bf4381 2355 //See header file for documentation
2356 AliGammaConversionAODObject omega;
77ac6f3e 2357 omega.SetPx(omegakf->GetPx());
2358 omega.SetPy(omegakf->GetPy());
2359 omega.SetPz(omegakf->GetPz());
04bf4381 2360 omega.SetChi2(omegakf->GetChi2());
77ac6f3e 2361 omega.SetE(omegakf->GetE());
04bf4381 2362 omega.SetIMass(mass);
2363 omega.SetLabel1(omegaDaughter);
6441e967 2364 // //dynamic_cast<AliAODPWG4Particle*>(fAODBranch->At(daughter1))->SetTagged(kTRUE);
04bf4381 2365 omega.SetLabel2(gammaDaughter);
6441e967 2366 AddToAODBranch(fAODOmega, omega);
04bf4381 2367}
2368
2369
2370
d7d7e825 2371void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
2372 // see header file for documentation
2373
6c84d371 2374 // for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){
2375 // for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){
037dc2db 2376
5ce758b0 2377 fESDEvent = fV0Reader->GetESDEvent();
2378
037dc2db 2379 if(fKFReconstructedGammasTClone->GetEntriesFast()>fV0Reader->GetNumberOfV0s()){
2380 cout<<"Warning, number of entries in the tclone is bigger than number of v0s"<<endl;
2381 }
2382
6c84d371 2383 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2384 for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();secondGammaIndex++){
d7d7e825 2385
6c84d371 2386 // AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];
2387 // AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];
a0b94e5c 2388
6c84d371 2389 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2390 AliKFParticle * twoGammaDecayCandidateDaughter1 = (AliKFParticle *)fKFReconstructedGammasTClone->At(secondGammaIndex);
a0b94e5c 2391
d7d7e825 2392 if(fElectronv1[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv1[firstGammaIndex]==fElectronv2[secondGammaIndex]){
2393 continue;
2394 }
2395 if(fElectronv2[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv2[firstGammaIndex]==fElectronv2[secondGammaIndex]){
2396 continue;
2397 }
a0b94e5c 2398
d7d7e825 2399 AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
2400
2401 Double_t massTwoGammaCandidate = 0.;
2402 Double_t widthTwoGammaCandidate = 0.;
2403 Double_t chi2TwoGammaCandidate =10000.;
2404 twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);
d707e3cf 2405 // if(twoGammaCandidate->GetNDF()>0){
2406 // chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();
2407 chi2TwoGammaCandidate = twoGammaCandidate->GetChi2();
d7d7e825 2408
70ef88b5 2409 fHistograms->FillHistogram("ESD_Mother_Chi2",chi2TwoGammaCandidate);
2e2da371 2410 if((chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
d7d7e825 2411
2412 TVector3 momentumVectorTwoGammaCandidate(twoGammaCandidate->GetPx(),twoGammaCandidate->GetPy(),twoGammaCandidate->GetPz());
2413 TVector3 spaceVectorTwoGammaCandidate(twoGammaCandidate->GetX(),twoGammaCandidate->GetY(),twoGammaCandidate->GetZ());
2414
2415 Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);
2416 Double_t rapidity;
96ade8ef 2417 if(twoGammaCandidate->GetE() - twoGammaCandidate->GetPz() <= 0 || twoGammaCandidate->GetE() + twoGammaCandidate->GetPz() <= 0){
2418 cout << "Error: |Pz| > E !!!! " << endl;
d7d7e825 2419 rapidity=0;
2420 }
2421 else{
2422 rapidity = 0.5*(TMath::Log((twoGammaCandidate->GetE() +twoGammaCandidate->GetPz()) / (twoGammaCandidate->GetE()-twoGammaCandidate->GetPz())));
2423 }
dc2883e4 2424
2425 if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut()){
2426 delete twoGammaCandidate;
2427 continue; // rapidity cut
2428 }
2429
2430
48682642 2431 Double_t alfa=0.0;
2432 if( (twoGammaDecayCandidateDaughter0->GetE()+twoGammaDecayCandidateDaughter1->GetE()) != 0){
2433 alfa=TMath::Abs((twoGammaDecayCandidateDaughter0->GetE()-twoGammaDecayCandidateDaughter1->GetE())
2434 /(twoGammaDecayCandidateDaughter0->GetE()+twoGammaDecayCandidateDaughter1->GetE()));
2435 }
d7d7e825 2436
1e7846f4 2437 if(openingAngleTwoGammaCandidate < fMinOpeningAngleGhostCut){
2438 delete twoGammaCandidate;
2439 continue; // minimum opening angle to avoid using ghosttracks
2440 }
2441
67381a40 2442 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
6de3471d 2443 fHistograms->FillHistogram("ESD_Mother_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);
2444 fHistograms->FillHistogram("ESD_Mother_Energy", twoGammaCandidate->GetE());
2445 fHistograms->FillHistogram("ESD_Mother_Pt", momentumVectorTwoGammaCandidate.Pt());
2446 fHistograms->FillHistogram("ESD_Mother_Eta", momentumVectorTwoGammaCandidate.Eta());
2447 fHistograms->FillHistogram("ESD_Mother_Rapid", rapidity);
2448 fHistograms->FillHistogram("ESD_Mother_Phi", spaceVectorTwoGammaCandidate.Phi());
2449 fHistograms->FillHistogram("ESD_Mother_Mass", massTwoGammaCandidate);
2450 fHistograms->FillHistogram("ESD_Mother_alfa", alfa);
9da64c46 2451 if(massTwoGammaCandidate>0.1 && massTwoGammaCandidate<0.15){
2452 fHistograms->FillHistogram("ESD_Mother_alfa_Pi0", alfa);
2453 }
dcdc851f 2454 if(massTwoGammaCandidate>0.5 && massTwoGammaCandidate<0.57){
2455 fHistograms->FillHistogram("ESD_Mother_alfa_Eta", alfa);
2456 }
2457
6de3471d 2458 fHistograms->FillHistogram("ESD_Mother_R", spaceVectorTwoGammaCandidate.Pt()); // Pt in Space == R!!!
2459 fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), spaceVectorTwoGammaCandidate.Pt());
2460 fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());
2461 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2462 fHistograms->FillHistogram("ESD_Mother_InvMass",massTwoGammaCandidate);
9c1cb6f7 2463 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2464 }
3c45d101 2465 if(alfa<0.1){
2466 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_E_alpha",massTwoGammaCandidate ,twoGammaCandidate->GetE());
2467 }
6de3471d 2468
9640a3d1 2469
9c1cb6f7 2470 if(fCalculateBackground){
5ce758b0 2471 /* Kenneth, just for testing*/
2472 AliGammaConversionBGHandler * bgHandlerTest = fV0Reader->GetBGHandler();
2473
2474 Int_t zbin= bgHandlerTest->GetZBinIndex(fV0Reader->GetVertexZ());
2475 Int_t mbin=0;
2476 Int_t multKAA=0;
2477 if(fUseTrackMultiplicityForBG == kTRUE){
2478 multKAA=fV0Reader->CountESDTracks();
2479 mbin = bgHandlerTest->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2480 }
2481 else{// means we use #v0s for multiplicity
2482 multKAA=fV0Reader->GetNGoodV0s();
2483 mbin = bgHandlerTest->GetMultiplicityBinIndex(fV0Reader->GetNGoodV0s());
2484 }
2485 // cout<<"Filling bin number "<<zbin<<" and "<<mbin<<endl;
2486 // cout<<"Corresponding to z = "<<fV0Reader->GetVertexZ()<<" and m = "<<multKAA<<endl;
6de3471d 2487 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2488 fHistograms->FillHistogram(Form("%d%dESD_Mother_InvMass",zbin,mbin),massTwoGammaCandidate);
2489 fHistograms->FillHistogram(Form("%d%dESD_Mother_InvMass_vs_Pt",zbin,mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2490 /* end Kenneth, just for testing*/
2491 fHistograms->FillHistogram(Form("%dESD_Mother_InvMass_vs_Pt",mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2492 }
5ce758b0 2493 }
2494 /* if(fCalculateBackground){
9c1cb6f7 2495 AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
2496 Int_t mbin= bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2497 fHistograms->FillHistogram(Form("%dESD_Mother_InvMass_vs_Pt",mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
5ce758b0 2498 }*/
ebcfaa7e 2499 // if(fDoNeutralMesonV0MCCheck){
2500 if(fDoMCTruth){
037dc2db 2501 //Kenneth: Checking the eta of the gamma to check the difference between 0.9 and 1.2
2502 Int_t indexKF1 = fKFReconstructedGammasV0Index.at(firstGammaIndex);
2503 if(indexKF1<fV0Reader->GetNumberOfV0s()){
2504 fV0Reader->GetV0(indexKF1);//updates to the correct v0
2505 Double_t eta1 = fV0Reader->GetMotherCandidateEta();
2506 Bool_t isRealPi0=kFALSE;
10e3319b 2507 Bool_t isRealEta=kFALSE;
037dc2db 2508 Int_t gamma1MotherLabel=-1;
2509 if(fV0Reader->HasSameMCMother() == kTRUE){
2510 //cout<<"This v0 is a real v0!!!!"<<endl;
2511 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
2512 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
2513 if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
2514 if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){
2515 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
2516 gamma1MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
2517 }
2518 }
23340c07 2519 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() ==111){
2520 gamma1MotherLabel=-111;
2521 }
d440c756 2522 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() ==221){
2523 gamma1MotherLabel=-221;
2524 }
037dc2db 2525 }
2526 }
2527 Int_t indexKF2 = fKFReconstructedGammasV0Index.at(secondGammaIndex);
2528 if(indexKF1 == indexKF2){
2529 cout<<"index of the two KF particles are the same.... should not happen"<<endl;
2530 }
2531 if(indexKF2<fV0Reader->GetNumberOfV0s()){
2532 fV0Reader->GetV0(indexKF2);
2533 Double_t eta2 = fV0Reader->GetMotherCandidateEta();
2534 Int_t gamma2MotherLabel=-1;
2535 if(fV0Reader->HasSameMCMother() == kTRUE){
2536 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
2537 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
2538 if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
2539 if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){
2540 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
2541 gamma2MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
2542 }
2543 }
23340c07 2544 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() ==111){
2545 gamma2MotherLabel=-111;
2546 }
d440c756 2547 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() ==221){
2548 gamma2MotherLabel=-221;
2549 }
2550
23340c07 2551 }
037dc2db 2552 }
2553 if(gamma1MotherLabel>=0 && gamma1MotherLabel==gamma2MotherLabel){
2554 if(fV0Reader->CheckIfPi0IsMother(gamma1MotherLabel)){
2555 isRealPi0=kTRUE;
2556 }
10e3319b 2557 if(fV0Reader->CheckIfEtaIsMother(gamma1MotherLabel)){
2558 isRealEta=kTRUE;
2559 }
2560
037dc2db 2561 }
6de3471d 2562 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2563 if(TMath::Abs(eta1)>0.9 && TMath::Abs(eta2)>0.9){
2564 // fHistograms->FillHistogram("ESD_Mother_InvMass_1212",massTwoGammaCandidate);
2565 // fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2566 if(isRealPi0 || isRealEta){
2567 fHistograms->FillHistogram("ESD_TruePi0_InvMass_1212",massTwoGammaCandidate);
2568 fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_1212",openingAngleTwoGammaCandidate);
2569 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2570 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2571 fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
d707e3cf 2572 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2573 }
6de3471d 2574
2575 if(!isRealPi0 && !isRealEta){
2576 if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
86e91247 2577 fHistograms->FillHistogram("ESD_TrueBckGG_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
6de3471d 2578 }else{
86e91247 2579 fHistograms->FillHistogram("ESD_TrueBckCont_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2580 }
d440c756 2581 if(gamma1MotherLabel==-111 || gamma2MotherLabel==-111 || gamma1MotherLabel==-221 || gamma2MotherLabel==-221){
86e91247 2582 fHistograms->FillHistogram("ESD_TruePi0DalitzCont_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
6de3471d 2583 }
bd6d9fa3 2584 }
bd6d9fa3 2585
6de3471d 2586 }
2587 else if(TMath::Abs(eta1)>0.9 || TMath::Abs(eta2)>0.9){
ebcfaa7e 2588 // fHistograms->FillHistogram("ESD_Mother_InvMass_0912",massTwoGammaCandidate);
2589 // fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
6de3471d 2590
2591 if(isRealPi0 || isRealEta){
2592 fHistograms->FillHistogram("ESD_TruePi0_InvMass_0912",massTwoGammaCandidate);
2593 fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_0912",openingAngleTwoGammaCandidate);
2594 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2595 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2596 fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
d707e3cf 2597 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2598 }
6de3471d 2599 if(!isRealPi0 && !isRealEta){
2600 if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
86e91247 2601 fHistograms->FillHistogram("ESD_TrueBckGG_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
6de3471d 2602 }else{
86e91247 2603 fHistograms->FillHistogram("ESD_TrueBckCont_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2604 }
d440c756 2605 if(gamma1MotherLabel==-111 || gamma2MotherLabel==-111 || gamma1MotherLabel==-221 || gamma2MotherLabel==-221){
86e91247 2606 fHistograms->FillHistogram("ESD_TruePi0DalitzCont_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
6de3471d 2607 }
bd6d9fa3 2608 }
2609 }
6de3471d 2610 else{
ebcfaa7e 2611 // fHistograms->FillHistogram("ESD_Mother_InvMass_0909",massTwoGammaCandidate);
2612 // fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
6de3471d 2613 if(isRealPi0 || isRealEta){
2614 fHistograms->FillHistogram("ESD_TruePi0_InvMass_0909",massTwoGammaCandidate);
2615 fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_0909",openingAngleTwoGammaCandidate);
2616 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2617 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2618 fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
d707e3cf 2619 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
6de3471d 2620 if(gamma1MotherLabel > fV0Reader->GetMCStack()->GetNprimary()){
86e91247 2621 fHistograms->FillHistogram("ESD_TruePi0Sec_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
6de3471d 2622 }
d707e3cf 2623 }
6de3471d 2624 if(!isRealPi0 && !isRealEta){
2625 if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
86e91247 2626 fHistograms->FillHistogram("ESD_TrueBckGG_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
6de3471d 2627 }else{
86e91247 2628 fHistograms->FillHistogram("ESD_TrueBckCont_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2629 }
d440c756 2630 if(gamma1MotherLabel==-111 || gamma2MotherLabel==-111 || gamma1MotherLabel==-221 || gamma2MotherLabel==-221 ){
86e91247 2631 fHistograms->FillHistogram("ESD_TruePi0DalitzCont_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
6de3471d 2632 }
bd6d9fa3 2633 }
037dc2db 2634 }
2635 }
2636 }
2637 }
2638 }
6de3471d 2639 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2640 if ( TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())<0.9 && TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())<0.9 ){
2641 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_Fiducial",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2642 fHistograms->FillHistogram("ESD_Mother_InvMass_Fiducial",massTwoGammaCandidate);
2643 }
2644
2645 if(TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())>0.9 && TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())>0.9){
2646 fHistograms->FillHistogram("ESD_Mother_InvMass_1212",massTwoGammaCandidate);
2647 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2648 }
2649 else if(TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())>0.9 || TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())>0.9){
2650 fHistograms->FillHistogram("ESD_Mother_InvMass_0912",massTwoGammaCandidate);
2651 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2652 }
2653 else{
2654 fHistograms->FillHistogram("ESD_Mother_InvMass_0909",massTwoGammaCandidate);
2655 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2656 }
ebcfaa7e 2657
6de3471d 2658 Double_t lowMassPi0=0.1;
2659 Double_t highMassPi0=0.15;
95203575 2660 if ( ( massTwoGammaCandidate > lowMassPi0) && (massTwoGammaCandidate < highMassPi0) ){
6de3471d 2661 new((*fKFReconstructedPi0sTClone)[fKFReconstructedPi0sTClone->GetEntriesFast()]) AliKFParticle(*twoGammaCandidate);
2662 fGammav1.push_back(firstGammaIndex);
2663 fGammav2.push_back(secondGammaIndex);
95203575 2664 if( fKFCreateAOD ) {
2665 AddPionToAOD(twoGammaCandidate, massTwoGammaCandidate, firstGammaIndex, secondGammaIndex);
2666 }
6de3471d 2667 }
ebcfaa7e 2668 }
6272370b 2669
2e2da371 2670 }
6de3471d 2671 delete twoGammaCandidate;
d7d7e825 2672 }
2673 }
2674}
2675
22eae5f9 2676void AliAnalysisTaskGammaConversion::AddPionToAOD(AliKFParticle * pionkf, Double_t mass, Int_t daughter1, Int_t daughter2) {
04bf4381 2677 //See header file for documentation
22eae5f9 2678 if(fOutputAODClassName.Contains("AODObject")) {
2679 AliGammaConversionAODObject pion;
2680 pion.SetPx(pionkf->GetPx());
2681 pion.SetPy(pionkf->GetPy());
2682 pion.SetPz(pionkf->GetPz());
2683 pion.SetChi2(pionkf->GetChi2());
2684 pion.SetE(pionkf->GetE());
2685 pion.SetIMass(mass);
2686 pion.SetLabel1(daughter1);
2687 pion.SetLabel2(daughter2);
2688 AddToAODBranch(fAODPi0, pion);
2689 } else {
2690 TLorentzVector momentum(pionkf->Px(),pionkf->Py(),pionkf->Pz(), pionkf->E());
2691 AliAODConversionParticle pion = AliAODConversionParticle(momentum);
2692 pion.SetTrackLabels( daughter1, daughter2 );
2693 pion.SetChi2(pionkf->GetChi2());
2694 AddToAODBranch(fAODPi0, pion);
2695
2696 }
04bf4381 2697}
2698
2699
48682642 2700/*
77ac6f3e 2701 void AliAnalysisTaskGammaConversion::ProcessConvPHOSGammasForNeutralMesonAnalysis(){
2702
6272370b 2703 // see header file for documentation
2704 // Analyse Pi0 with one photon from Phos and 1 photon from conversions
2705
2706
2707
2708 Double_t vtx[3];
2709 vtx[0] = fV0Reader->GetPrimaryVertex()->GetX();
2710 vtx[1] = fV0Reader->GetPrimaryVertex()->GetY();
2711 vtx[2] = fV0Reader->GetPrimaryVertex()->GetZ();
2712
2713
2714 // Loop over all CaloClusters and consider only the PHOS ones:
2715 AliESDCaloCluster *clu;
2716 TLorentzVector pPHOS;
2717 TLorentzVector gammaPHOS;
2718 TLorentzVector gammaGammaConv;
2719 TLorentzVector pi0GammaConvPHOS;
2720 TLorentzVector gammaGammaConvBck;
2721 TLorentzVector pi0GammaConvPHOSBck;
2722
2723
2724 for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
77ac6f3e 2725 clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
2726 if ( !clu->IsPHOS() || clu->E()<0.1 ) continue;
2727 clu ->GetMomentum(pPHOS ,vtx);
2728 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2729 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2730 gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz(),0.);
2731 gammaPHOS.SetXYZM(pPHOS.Px(),pPHOS.Py(),pPHOS.Pz(),0.);
2732 pi0GammaConvPHOS=gammaGammaConv+gammaPHOS;
2733 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS",pi0GammaConvPHOS.M());
2734 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvPHOS",pi0GammaConvPHOS.M(),pi0GammaConvPHOS.Pt());
6272370b 2735
77ac6f3e 2736 TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz());
2737 TVector3 v3D1(gammaPHOS.Px(),gammaPHOS.Py(),gammaPHOS.Pz());
2738 Double_t opanConvPHOS= v3D0.Angle(v3D1);
2739 if ( opanConvPHOS < 0.35){
2740 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanLow",pi0GammaConvPHOS.M());
2741 }else{
2742 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanHigh",pi0GammaConvPHOS.M());
2743 }
6272370b 2744
77ac6f3e 2745 }
2746
2747 // Now the LorentVector pPHOS is obtained and can be paired with the converted proton
6272370b 2748 }
2749 //==== End of the PHOS cluster selection ============
2750 TLorentzVector pEMCAL;
2751 TLorentzVector gammaEMCAL;
2752 TLorentzVector pi0GammaConvEMCAL;
2753 TLorentzVector pi0GammaConvEMCALBck;
2754
77ac6f3e 2755 for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
2756 clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
2757 if ( !clu->IsEMCAL() || clu->E()<0.1 ) continue;
2758 if (clu->GetNCells() <= 1) continue;
2759 if ( clu->GetTOF()*1e9 < 550 || clu->GetTOF()*1e9 > 750) continue;
6272370b 2760
77ac6f3e 2761 clu ->GetMomentum(pEMCAL ,vtx);
2762 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2763 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2764 gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),
2765 twoGammaDecayCandidateDaughter0->Py(),
2766 twoGammaDecayCandidateDaughter0->Pz(),0.);
2767 gammaEMCAL.SetXYZM(pEMCAL.Px(),pEMCAL.Py(),pEMCAL.Pz(),0.);
2768 pi0GammaConvEMCAL=gammaGammaConv+gammaEMCAL;
2769 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL",pi0GammaConvEMCAL.M());
2770 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL",pi0GammaConvEMCAL.M(),pi0GammaConvEMCAL.Pt());
2771 TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),
2772 twoGammaDecayCandidateDaughter0->Py(),
2773 twoGammaDecayCandidateDaughter0->Pz());
2774 TVector3 v3D1(gammaEMCAL.Px(),gammaEMCAL.Py(),gammaEMCAL.Pz());
6272370b 2775
77ac6f3e 2776
2777 Double_t opanConvEMCAL= v3D0.Angle(v3D1);
2778 if ( opanConvEMCAL < 0.35){
2779 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanLow",pi0GammaConvEMCAL.M());
2780 }else{
2781 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanHigh",pi0GammaConvEMCAL.M());
2782 }
2783
2784 }
2785 if(fCalculateBackground){
2786 for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2787 AliGammaConversionKFVector * previousEventV0s = fV0Reader->GetBGGoodV0s(nEventsInBG);
2788 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2789 AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2790 gammaGammaConvBck.SetXYZM(previousGoodV0.Px(),
2791 previousGoodV0.Py(),
2792 previousGoodV0.Pz(),0.);
2793 pi0GammaConvEMCALBck=gammaGammaConvBck+gammaEMCAL;
2794 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M());
2795 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M(),
2796 pi0GammaConvEMCALBck.Pt());
2797 }
2798 }
48682642 2799
77ac6f3e 2800 // Now the LorentVector pEMCAL is obtained and can be paired with the converted proton
2801 } // end of checking if background photons are available
2802 }
6272370b 2803 //==== End of the PHOS cluster selection ============
77ac6f3e 2804
2805 }
48682642 2806*/
6272370b 2807
77ac6f3e 2808void AliAnalysisTaskGammaConversion::MoveParticleAccordingToVertex(AliKFParticle * particle,const AliGammaConversionBGHandler::GammaConversionVertex *vertex){
5ce758b0 2809 //see header file for documentation
2810
77ac6f3e 2811 Double_t dx = vertex->fX - fESDEvent->GetPrimaryVertex()->GetX();
2812 Double_t dy = vertex->fY - fESDEvent->GetPrimaryVertex()->GetY();
2813 Double_t dz = vertex->fZ - fESDEvent->GetPrimaryVertex()->GetZ();
5ce758b0 2814
2815 // cout<<"dx, dy, dz: ["<<dx<<","<<dy<<","<<dz<<"]"<<endl;
77ac6f3e 2816 particle->X() = particle->GetX() - dx;
2817 particle->Y() = particle->GetY() - dy;
2818 particle->Z() = particle->GetZ() - dz;
5ce758b0 2819}
2820
111d75df 2821void AliAnalysisTaskGammaConversion::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle){
5510ebe1 2822 // Before rotate needs to be moved to position 0,0,0, ; move back after rotation
2823 Double_t dx = fESDEvent->GetPrimaryVertex()->GetX()-0.;
2824 Double_t dy = fESDEvent->GetPrimaryVertex()->GetY()-0.;
2825 Double_t dz = fESDEvent->GetPrimaryVertex()->GetZ()-0.;
2826
2827 kfParticle->X() = kfParticle->GetX() - dx;
2828 kfParticle->Y() = kfParticle->GetY() - dy;
2829 kfParticle->Z() = kfParticle->GetZ() - dz;
2830
2831
77ac6f3e 2832 // Rotate the kf particle
111d75df 2833 Double_t c = cos(angle);
2834 Double_t s = sin(angle);
2835
5510ebe1 2836 Double_t mA[8][ 8];
2837 for( Int_t i=0; i<8; i++ ){
2838 for( Int_t j=0; j<8; j++){
77ac6f3e 2839 mA[i][j] = 0;
111d75df 2840 }
2841 }
5510ebe1 2842 for( int i=0; i<8; i++ ){
77ac6f3e 2843 mA[i][i] = 1;
111d75df 2844 }
77ac6f3e 2845 mA[0][0] = c; mA[0][1] = s;
2846 mA[1][0] = -s; mA[1][1] = c;
2847 mA[3][3] = c; mA[3][4] = s;
2848 mA[4][3] = -s; mA[4][4] = c;
111d75df 2849
5510ebe1 2850 Double_t mAC[8][8];
2851 Double_t mAp[8];
111d75df 2852
5510ebe1 2853 for( Int_t i=0; i<8; i++ ){
77ac6f3e 2854 mAp[i] = 0;
5510ebe1 2855 for( Int_t k=0; k<8; k++){
77ac6f3e 2856 mAp[i]+=mA[i][k] * kfParticle->GetParameter(k);
111d75df 2857 }
2858 }
2859
5510ebe1 2860 for( Int_t i=0; i<8; i++){
77ac6f3e 2861 kfParticle->Parameter(i) = mAp[i];
111d75df 2862 }
2863
5510ebe1 2864 for( Int_t i=0; i<8; i++ ){
2865 for( Int_t j=0; j<8; j++ ){
77ac6f3e 2866 mAC[i][j] = 0;
5510ebe1 2867 for( Int_t k=0; k<8; k++ ){
77ac6f3e 2868 mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j);
111d75df 2869 }
2870 }
2871 }
2872
5510ebe1 2873 for( Int_t i=0; i<8; i++ ){
111d75df 2874 for( Int_t j=0; j<=i; j++ ){
2875 Double_t xx = 0;
5510ebe1 2876 for( Int_t k=0; k<8; k++){
77ac6f3e 2877 xx+= mAC[i][k]*mA[j][k];
111d75df 2878 }
2879 kfParticle->Covariance(i,j) = xx;
2880 }
2881 }
5510ebe1 2882
2883 Double_t dx1 = 0.-fESDEvent->GetPrimaryVertex()->GetX();
2884 Double_t dy1 = 0.-fESDEvent->GetPrimaryVertex()->GetY();
2885 Double_t dz1 = 0.-fESDEvent->GetPrimaryVertex()->GetZ();
2886
2887 kfParticle->X() = kfParticle->GetX() - dx1;
2888 kfParticle->Y() = kfParticle->GetY() - dy1;
2889 kfParticle->Z() = kfParticle->GetZ() - dz1;
2890
111d75df 2891}
2892
2893
d7d7e825 2894void AliAnalysisTaskGammaConversion::CalculateBackground(){
2895 // see header file for documentation
5e55d806 2896
2897
2898 TClonesArray * currentEventV0s = fV0Reader->GetCurrentEventGoodV0s();
2899
5ce758b0 2900 AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
2901
2902 Int_t zbin= bgHandler->GetZBinIndex(fV0Reader->GetVertexZ());
2903 Int_t mbin = 0;
2904 if(fUseTrackMultiplicityForBG == kTRUE){
2905 mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2906 }
2907 else{
2908 mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->GetNGoodV0s());
2909 }
2910
111d75df 2911 if(fDoRotation == kTRUE){
5510ebe1 2912 TRandom3 *random = new TRandom3(0);
2913
111d75df 2914 for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2915 AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
2916 for(Int_t iCurrent2=iCurrent+1;iCurrent2<currentEventV0s->GetEntriesFast();iCurrent2++){
77ac6f3e 2917 for(Int_t nRandom=0;nRandom<fNRandomEventsForBG;nRandom++){
111d75df 2918
2919 AliKFParticle currentEventGoodV02 = *(AliKFParticle *)(currentEventV0s->At(iCurrent2));
3bfbe89a 2920
2921 if(fCheckBGProbability == kTRUE){
2922 Double_t massBGprob =0.;
2923 Double_t widthBGprob = 0.;
2924 AliKFParticle *backgroundCandidateProb = new AliKFParticle(currentEventGoodV0,currentEventGoodV02);
2925 backgroundCandidateProb->GetMass(massBGprob,widthBGprob);
2926 if(massBGprob>0.1 && massBGprob<0.14){
2927 if(random->Rndm()>bgHandler->GetBGProb(zbin,mbin)){
2928 delete backgroundCandidateProb;
2929 continue;
2930 }
2931 }
2932 delete backgroundCandidateProb;
2933 }
111d75df 2934
77ac6f3e 2935 Double_t nRadiansPM = fNDegreesPMBackground*TMath::Pi()/180;
111d75df 2936
2937 Double_t rotationValue = random->Rndm()*2*nRadiansPM + TMath::Pi()-nRadiansPM;
5510ebe1 2938
111d75df 2939 RotateKFParticle(&currentEventGoodV02,rotationValue);
5ce758b0 2940
111d75df 2941 AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,currentEventGoodV02);
5e55d806 2942
5ce758b0 2943 Double_t massBG =0.;
2944 Double_t widthBG = 0.;
2945 Double_t chi2BG =10000.;
2946 backgroundCandidate->GetMass(massBG,widthBG);
111d75df 2947 // if(backgroundCandidate->GetNDF()>0){
2948 chi2BG = backgroundCandidate->GetChi2();
2949 if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
2950
2951 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2952 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2953
2954 Double_t openingAngleBG = currentEventGoodV0.GetAngle(currentEventGoodV02);
2955
2956 Double_t rapidity;
2957 if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
2958 else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2959
dc2883e4 2960 if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ){
2961 delete backgroundCandidate;
2962 continue; // rapidity cut
2963 }
2964
111d75df 2965
2966 Double_t alfa=0.0;
2967 if( (currentEventGoodV0.GetE()+currentEventGoodV02.GetE()) != 0){
2968 alfa=TMath::Abs((currentEventGoodV0.GetE()-currentEventGoodV02.GetE())
2969 /(currentEventGoodV0.GetE()+currentEventGoodV02.GetE()));
2970 }
2971
2972
2973 if(openingAngleBG < fMinOpeningAngleGhostCut ){
2974 delete backgroundCandidate;
2975 continue; // minimum opening angle to avoid using ghosttracks
2976 }
2977
2978 // original
111d75df 2979 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
6de3471d 2980 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2981 fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2982 fHistograms->FillHistogram("ESD_Background_Pt", momentumVectorbackgroundCandidate.Pt());
2983 fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2984 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2985 fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2986 fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2987 fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2988 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2989 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2990 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2991 fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
111d75df 2992 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
6de3471d 2993
dcdc851f 2994 if(massBG>0.1 && massBG<0.15){
2995 fHistograms->FillHistogram("ESD_Background_alfa_Pi0", alfa);
2996 }
2997 if(massBG>0.5 && massBG<0.57){
2998 fHistograms->FillHistogram("ESD_Background_alfa_Eta", alfa);
2999 }
6de3471d 3000
3001 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(currentEventGoodV02.GetEta())<0.9 ){
3002 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
3003 fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
3004 }
3005
3006 fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
3007 fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
3008 fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin), momentumVectorbackgroundCandidate.Pt());
3009 fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
3010 fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
3011 fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
3012 fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
3013 fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
3014 fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
3015 fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
3016 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
3017 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
3018
3019 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(currentEventGoodV02.GetEta())<0.9 ){
3020 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
3021 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
3022 }
111d75df 3023 }
3c45d101 3024 if(alfa<0.1){
3025 fHistograms->FillHistogram("ESD_Background_InvMass_vs_E_alpha",massBG ,backgroundCandidate->GetE());
3026 }
3027
2e2da371 3028 }
111d75df 3029 //}
5ce758b0 3030 delete backgroundCandidate;
3031 }
3032 }
3033 }
111d75df 3034 delete random;
5ce758b0 3035 }
111d75df 3036 else{ // means no rotation
3037 AliGammaConversionBGHandler::GammaConversionVertex *bgEventVertex = NULL;
3038
3039 if(fUseTrackMultiplicityForBG){
3040 // cout<<"Using charged track multiplicity for background calculation"<<endl;
3041 for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
5ce758b0 3042
111d75df 3043 AliGammaConversionKFVector * previousEventV0s = bgHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);//fV0Reader->GetBGGoodV0s(nEventsInBG);
3044
5ce758b0 3045 if(fMoveParticleAccordingToVertex == kTRUE){
3046 bgEventVertex = bgHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
3047 }
3048
3049 for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
3050 AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
3051 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
3052 AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
111d75df 3053 AliKFParticle previousGoodV0test = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
5ce758b0 3054
111d75df 3055 //cout<<"Primary Vertex event: ["<<fESDEvent->GetPrimaryVertex()->GetX()<<","<<fESDEvent->GetPrimaryVertex()->GetY()<<","<<fESDEvent->GetPrimaryVertex()->GetZ()<<"]"<<endl;
3056 //cout<<"BG prim Vertex event: ["<<bgEventVertex->fX<<","<<bgEventVertex->fY<<","<<bgEventVertex->fZ<<"]"<<endl;
3057
3058 //cout<<"XYZ of particle before transport: ["<<previousGoodV0.X()<<","<<previousGoodV0.Y()<<","<<previousGoodV0.Z()<<"]"<<endl;
5ce758b0 3059 if(fMoveParticleAccordingToVertex == kTRUE){
3060 MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
3061 }
111d75df 3062 //cout<<"XYZ of particle after transport: ["<<previousGoodV0.X()<<","<<previousGoodV0.Y()<<","<<previousGoodV0.Z()<<"]"<<endl;
5ce758b0 3063
3064 AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
111d75df 3065
5ce758b0 3066 Double_t massBG =0.;
3067 Double_t widthBG = 0.;
3068 Double_t chi2BG =10000.;
3069 backgroundCandidate->GetMass(massBG,widthBG);
5510ebe1 3070
111d75df 3071 // if(backgroundCandidate->GetNDF()>0){
3072 // chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
3073 chi2BG = backgroundCandidate->GetChi2();
3074 if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
3075
3076 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
3077 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
3078
3079 Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
3080
3081 Double_t rapidity;
3082
3083 if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() <= 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() <= 0){
3084 cout << "Error: |Pz| > E !!!! " << endl;
3085 rapidity=0;
3086 } else {
3087 rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
3088 }
dc2883e4 3089 if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ){
3090 delete backgroundCandidate;
3091 continue; // rapidity cut
3092 }
3093
3094
111d75df 3095 Double_t alfa=0.0;
3096 if( (currentEventGoodV0.GetE()+previousGoodV0.GetE()) != 0){
3097 alfa=TMath::Abs((currentEventGoodV0.GetE()-previousGoodV0.GetE())
3098 /(currentEventGoodV0.GetE()+previousGoodV0.GetE()));
3099 }
3100
3101
3102 if(openingAngleBG < fMinOpeningAngleGhostCut ){
3103 delete backgroundCandidate;
3104 continue; // minimum opening angle to avoid using ghosttracks
3105 }
3106
3107 // original
111d75df 3108 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
6de3471d 3109 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
3110 fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
3111 fHistograms->FillHistogram("ESD_Background_Pt", momentumVectorbackgroundCandidate.Pt());
3112 fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
3113 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
3114 fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
3115 fHistograms->FillHistogram("ESD_Background_Mass", massBG);
3116 fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
3117 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
3118 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
3119 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
3120 fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
111d75df 3121 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
6de3471d 3122
dcdc851f 3123 if(massBG>0.1 && massBG<0.15){
3124 fHistograms->FillHistogram("ESD_Background_alfa_Pi0", alfa);
3125 }
3126 if(massBG>0.5 && massBG<0.57){
3127 fHistograms->FillHistogram("ESD_Background_alfa_Eta", alfa);
3128 }
6de3471d 3129
3130 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
3131 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
3132 fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
3133 }
3134
3135 // test
3136 fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
3137 fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
3138 fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin), momentumVectorbackgroundCandidate.Pt());
3139 fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
3140 fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
3141 fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
3142 fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
3143 fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
3144 fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
3145 fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
3146 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
3147 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
3148
3149 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
3150 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
3151 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
3152 }
3153 // }
111d75df 3154 }
3c45d101 3155 if(alfa<0.1){
3156 fHistograms->FillHistogram("ESD_Background_InvMass_vs_E_alpha",massBG ,backgroundCandidate->GetE());
3157 }
3158
111d75df 3159 }
3160 delete backgroundCandidate;
3161 }
3162 }
3163 }
3164 }
3165 else{ // means using #V0s for multiplicity
3166
3167 // cout<<"Using the v0 multiplicity to calculate background"<<endl;
3168
3169 fHistograms->FillHistogram("ESD_Background_z_m",zbin,mbin);
3170 fHistograms->FillHistogram("ESD_Mother_multpilicityVSv0s",fV0Reader->CountESDTracks(),fV0Reader->GetNumberOfV0s());
3171
3172 for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
3173 AliGammaConversionKFVector * previousEventV0s = bgHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);// fV0Reader->GetBGGoodV0s(nEventsInBG);
3174 if(previousEventV0s){
3175
3176 if(fMoveParticleAccordingToVertex == kTRUE){
3177 bgEventVertex = bgHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
3178 }
3179
3180 for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
3181 AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
3182 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
3183 AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
3184
3185 if(fMoveParticleAccordingToVertex == kTRUE){
3186 MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
5ce758b0 3187 }
111d75df 3188
3189 AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
3190 Double_t massBG =0.;
3191 Double_t widthBG = 0.;
3192 Double_t chi2BG =10000.;
3193 backgroundCandidate->GetMass(massBG,widthBG);
5510ebe1 3194
111d75df 3195 /* if(backgroundCandidate->GetNDF()>0){
3196 chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
3197 {//remember to remove
3198 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
3199 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
3200
3201 Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
3202 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle_nochi2", openingAngleBG);
3203 }
3204 */
2e2da371 3205 chi2BG = backgroundCandidate->GetChi2();
3206 if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
5ce758b0 3207 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
3208 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
3209
3210 Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
3211
3212 Double_t rapidity;
3213 if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
3214 else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
3215
dc2883e4 3216 if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ){
3217 delete backgroundCandidate;
3218 continue; // rapidity cut
3219 }
3220
3221
5ce758b0 3222 Double_t alfa=0.0;
3223 if( (currentEventGoodV0.GetE()+previousGoodV0.GetE()) != 0){
3224 alfa=TMath::Abs((currentEventGoodV0.GetE()-previousGoodV0.GetE())
3225 /(currentEventGoodV0.GetE()+previousGoodV0.GetE()));
3226 }
3227
3228
3229 if(openingAngleBG < fMinOpeningAngleGhostCut ){
3230 delete backgroundCandidate;
3231 continue; // minimum opening angle to avoid using ghosttracks
3232 }
3233
5ce758b0 3234 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
6de3471d 3235 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
3236 fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
3237 fHistograms->FillHistogram("ESD_Background_Pt", momentumVectorbackgroundCandidate.Pt());
3238 fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
3239 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
3240 fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
3241 fHistograms->FillHistogram("ESD_Background_Mass", massBG);
3242 fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
3243 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
3244 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
3245 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
3246 fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
3247
3248
5ce758b0 3249 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
3c45d101 3250
dcdc851f 3251 if(massBG>0.1 && massBG<0.15){
3252 fHistograms->FillHistogram("ESD_Background_alfa_Pi0", alfa);
3253 }
3254 if(massBG>0.5 && massBG<0.57){
3255 fHistograms->FillHistogram("ESD_Background_alfa_Eta", alfa);
3256 }
3257
6de3471d 3258 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
3259 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
3260 fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
3261 }
3262
3263 if(massBG>0.5 && massBG<0.6){
3264 fHistograms->FillHistogram("ESD_Background_alfa_pt0506",momentumVectorbackgroundCandidate.Pt(),alfa);
3265 }
3266 if(massBG>0.3 && massBG<0.4){
3267 fHistograms->FillHistogram("ESD_Background_alfa_pt0304",momentumVectorbackgroundCandidate.Pt(),alfa);
3268 }
3269
3270 // test
3271 fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
3272 fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
3273 fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin), momentumVectorbackgroundCandidate.Pt());
3274 fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
3275 fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
3276 fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
3277 fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
3278 fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
3279 fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
3280 fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
3281 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
3282 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
3283
3284 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
3285 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
3286 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
3287 }
5ce758b0 3288 }
5ce758b0 3289
6de3471d 3290 if(alfa<0.1){
3291 fHistograms->FillHistogram("ESD_Background_InvMass_vs_E_alpha",massBG ,backgroundCandidate->GetE());
5ce758b0 3292 }
2e2da371 3293 // }
5ce758b0 3294 }
111d75df 3295 delete backgroundCandidate;
3296 }
5ce758b0 3297 }
3298 }
d7d7e825 3299 }
111d75df 3300 } // end else (means use #v0s as multiplicity)
3301 } // end no rotation
d7d7e825 3302}
3303
3304
d7d7e825 3305void AliAnalysisTaskGammaConversion::ProcessGammasForGammaJetAnalysis(){
3306 //ProcessGammasForGammaJetAnalysis
a0b94e5c 3307
d7d7e825 3308 Double_t distIsoMin;
a0b94e5c 3309
d7d7e825 3310 CreateListOfChargedParticles();
a0b94e5c 3311
3312
6c84d371 3313 // for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
3314 for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
3315 AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
01b7fdcc 3316 TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
a0b94e5c 3317
01b7fdcc 3318 if( momentumVectorCurrentGamma.Pt()> fMinPtForGammaJet){
d7d7e825 3319 distIsoMin=GetMinimumDistanceToCharge(gammaIndex);
a0b94e5c 3320
d7d7e825 3321 if (distIsoMin > fMinIsoConeSize && fLeadingChargedIndex>=0){
01b7fdcc 3322 CalculateJetCone(gammaIndex);
d7d7e825 3323 }
3324 }
3325 }
3326}
3327
6272370b 3328//____________________________________________________________________
77ac6f3e 3329Bool_t AliAnalysisTaskGammaConversion::IsGoodImpPar(const AliESDtrack *const track)
6272370b 3330{
3331//
3332// check whether particle has good DCAr(Pt) impact
3333// parameter. Only for TPC+ITS tracks (7*sigma cut)
3334// Origin: Andrea Dainese
3335//
3336
3337Float_t d0z0[2],covd0z0[3];
3338track->GetImpactParameters(d0z0,covd0z0);
3339Float_t sigma= 0.0050+0.0060/TMath::Power(track->Pt(),0.9);
3340Float_t d0max = 7.*sigma;
3341if(TMath::Abs(d0z0[0]) < d0max) return kTRUE;
3342
3343return kFALSE;
3344}
3345
3346
d7d7e825 3347void AliAnalysisTaskGammaConversion::CreateListOfChargedParticles(){
3348 // CreateListOfChargedParticles
a0b94e5c 3349
d7d7e825 3350 fESDEvent = fV0Reader->GetESDEvent();
037dc2db 3351 Int_t numberOfESDTracks=0;
d7d7e825 3352 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
3353 AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
a0b94e5c 3354
d7d7e825 3355 if(!curTrack){
3356 continue;
3357 }
d707e3cf 3358 // Not needed if Standard function used.
3359// if(!IsGoodImpPar(curTrack)){
3360// continue;
3361// }
a0b94e5c 3362
d7d7e825 3363 if(fEsdTrackCuts->AcceptTrack(curTrack) ){
6c84d371 3364 new((*fChargedParticles)[fChargedParticles->GetEntriesFast()]) AliESDtrack(*curTrack);
3365 // fChargedParticles.push_back(curTrack);
d7d7e825 3366 fChargedParticlesId.push_back(iTracks);
037dc2db 3367 numberOfESDTracks++;
d7d7e825 3368 }
3369 }
6746e1e1 3370// Moved to UserExec using CountAcceptedTracks function. runjet is not needed by default
3371// fHistograms->FillHistogram("ESD_NumberOfGoodESDTracks",numberOfESDTracks);
3372// cout<<"esdtracks::"<< numberOfESDTracks<<endl;
3373// if (fV0Reader->GetNumberOfContributorsVtx()>=1){
3374// fHistograms->FillHistogram("ESD_NumberOfGoodESDTracksVtx",numberOfESDTracks);
3375// }
d7d7e825 3376}
9c1cb6f7 3377void AliAnalysisTaskGammaConversion::RecalculateV0ForGamma(){
77ac6f3e 3378 //recalculates v0 for gamma
3379
9c1cb6f7 3380 Double_t massE=0.00051099892;
3381 TLorentzVector curElecPos;
3382 TLorentzVector curElecNeg;
3383 TLorentzVector curGamma;
3384
3385 TLorentzVector curGammaAt;
3386 TLorentzVector curElecPosAt;
3387 TLorentzVector curElecNegAt;
3388 AliKFVertex primVtxGamma(*(fESDEvent->GetPrimaryVertex()));
3389 AliKFVertex primVtxImprovedGamma = primVtxGamma;
3390
3391 const AliESDVertex *vtxT3D=fESDEvent->GetPrimaryVertex();
3392
3393 Double_t xPrimaryVertex=vtxT3D->GetXv();
3394 Double_t yPrimaryVertex=vtxT3D->GetYv();
3395 Double_t zPrimaryVertex=vtxT3D->GetZv();
70ef88b5 3396 // Float_t primvertex[3]={xPrimaryVertex,yPrimaryVertex,zPrimaryVertex};
9c1cb6f7 3397
3398 Float_t nsigmaTPCtrackPos;
3399 Float_t nsigmaTPCtrackNeg;
3400 Float_t nsigmaTPCtrackPosToPion;
3401 Float_t nsigmaTPCtrackNegToPion;
3402 AliKFParticle* negKF=NULL;
3403 AliKFParticle* posKF=NULL;
3404
3405 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
3406 AliESDtrack* posTrack = fESDEvent->GetTrack(iTracks);
3407 if(!posTrack){
3408 continue;
3409 }
3410 if (posKF) delete posKF; posKF=NULL;
3411 if(posTrack->GetSign()<0) continue;
3412 if(!(posTrack->GetStatus() & AliESDtrack::kTPCrefit))continue;
3413 if(posTrack->GetKinkIndex(0)>0 ) continue;
3414 if(posTrack->GetNcls(1)<50)continue;
3415 Double_t momPos[3];
3416 // posTrack->GetConstrainedPxPyPz(momPos);
3417 posTrack->GetPxPyPz(momPos);
3418 AliESDtrack *ptrk=fESDEvent->GetTrack(iTracks);
3419 curElecPos.SetXYZM(momPos[0],momPos[1],momPos[2],massE);
3420 if(TMath::Abs(curElecPos.Eta())<0.9) continue;
3421 posKF = new AliKFParticle( *(posTrack),-11);
3422
3423 nsigmaTPCtrackPos = fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kElectron);
3424 nsigmaTPCtrackPosToPion = fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kPion);
3425
3426 if ( nsigmaTPCtrackPos>5.|| nsigmaTPCtrackPos<-2.){
3427 continue;
3428 }
3429
3430 if(pow((momPos[0]*momPos[0]+momPos[1]*momPos[1]+momPos[2]*momPos[2]),0.5)>0.5 && nsigmaTPCtrackPosToPion<1){
3431 continue;
3432 }
3433
3434
3435
3436 for(Int_t jTracks = 0; jTracks < fESDEvent->GetNumberOfTracks(); jTracks++){
3437 AliESDtrack* negTrack = fESDEvent->GetTrack(jTracks);
3438 if(!negTrack){
3439 continue;
3440 }
3441 if (negKF) delete negKF; negKF=NULL;
3442 if(negTrack->GetSign()>0) continue;
3443 if(!(negTrack->GetStatus() & AliESDtrack::kTPCrefit))continue;
3444 if(negTrack->GetKinkIndex(0)>0 ) continue;
3445 if(negTrack->GetNcls(1)<50)continue;
3446 Double_t momNeg[3];
3447 // negTrack->GetConstrainedPxPyPz(momNeg);
3448 negTrack->GetPxPyPz(momNeg);
3449
3450 nsigmaTPCtrackNeg = fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kElectron);
3451 nsigmaTPCtrackNegToPion = fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kPion);
3452 if ( nsigmaTPCtrackNeg>5. || nsigmaTPCtrackNeg<-2.){
3453 continue;
3454 }
3455 if(pow((momNeg[0]*momNeg[0]+momNeg[1]*momNeg[1]+momNeg[2]*momNeg[2]),0.5)>0.5 && nsigmaTPCtrackNegToPion<1){
3456 continue;
3457 }
3458 AliESDtrack *ntrk=fESDEvent->GetTrack(jTracks);
3459 curElecNeg.SetXYZM(momNeg[0],momNeg[1],momNeg[2],massE);
3460 if(TMath::Abs(curElecNeg.Eta())<0.9) continue;
3461 negKF = new AliKFParticle( *(negTrack) ,11);
3462
3463 Double_t b=fESDEvent->GetMagneticField();
3464 Double_t xn, xp, dca=ntrk->GetDCA(ptrk,b,xn,xp);
3465 AliExternalTrackParam nt(*ntrk), pt(*ptrk);
3466 nt.PropagateTo(xn,b); pt.PropagateTo(xp,b);
3467
3468
3469 //--- Like in ITSV0Finder
3470 AliExternalTrackParam ntAt0(*ntrk), ptAt0(*ptrk);
3471 Double_t xxP,yyP,alphaP;
3472 Double_t rP[3];
3473
3474 // if (!ptAt0.GetGlobalXYZat(ptAt0->GetX(),xxP,yyP,zzP)) continue;
3475 if (!ptAt0.GetXYZAt(ptAt0.GetX(),b,rP)) continue;
3476 xxP=rP[0];
3477 yyP=rP[1];
3478 alphaP = TMath::ATan2(yyP,xxP);
3479
3480
3481 ptAt0.Propagate(alphaP,0,b);
3482 Float_t ptfacP = (1.+100.*TMath::Abs(ptAt0.GetC(b)));
3483
3484 // Double_t distP = ptAt0.GetY();
3485 Double_t normP = ptfacP*TMath::Sqrt(ptAt0.GetSigmaY2());
3486 Double_t normdist0P = TMath::Abs(ptAt0.GetY()/normP);
3487 Double_t normdist1P = TMath::Abs((ptAt0.GetZ()-zPrimaryVertex)/(ptfacP*TMath::Sqrt(ptAt0.GetSigmaZ2())));
3488 Double_t normdistP = TMath::Sqrt(normdist0P*normdist0P+normdist1P*normdist1P);
3489
3490
3491 Double_t xxN,yyN,alphaN;
3492 Double_t rN[3];
3493 // if (!ntAt0.GetGlobalXYZat(ntAt0->GetX(),xxN,yyN,zzN)) continue;
3494 if (!ntAt0.GetXYZAt(ntAt0.GetX(),b,rN)) continue;
3495 xxN=rN[0];
3496 yyN=rN[1];
3497
3498 alphaN = TMath::ATan2(yyN,xxN);
3499
3500 ntAt0.Propagate(alphaN,0,b);
3501
3502 Float_t ptfacN = (1.+100.*TMath::Abs(ntAt0.GetC(b)));
3503 // Double_t distN = ntAt0.GetY();
3504 Double_t normN = ptfacN*TMath::Sqrt(ntAt0.GetSigmaY2());
3505 Double_t normdist0N = TMath::Abs(ntAt0.GetY()/normN);
3506 Double_t normdist1N = TMath::Abs((ntAt0.GetZ()-zPrimaryVertex)/(ptfacN*TMath::Sqrt(ntAt0.GetSigmaZ2())));
3507 Double_t normdistN = TMath::Sqrt(normdist0N*normdist0N+normdist1N*normdist1N);
3508
3509 //-----------------------------
3510
3511 Double_t momNegAt[3];
3512 nt.GetPxPyPz(momNegAt);
3513 curElecNegAt.SetXYZM(momNegAt[0],momNegAt[1],momNegAt[2],massE);
3514
3515 Double_t momPosAt[3];
3516 pt.GetPxPyPz(momPosAt);
3517 curElecPosAt.SetXYZM(momPosAt[0],momPosAt[1],momPosAt[2],massE);
3518 if(dca>1){
3519 continue;
3520 }
3521
3522 // Double_t dneg= negTrack->GetD(xPrimaryVertex,yPrimaryVertex,b);
3523 // Double_t dpos= posTrack->GetD(xPrimaryVertex,yPrimaryVertex,b);
3524 AliESDv0 vertex(nt,jTracks,pt,iTracks);
3525
3526
3527 Float_t cpa=vertex.GetV0CosineOfPointingAngle(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex);
3528
9c1cb6f7 3529
70ef88b5 3530
9c1cb6f7 3531 // cout<< "v0Rr::"<< v0Rr<<endl;
3532 // if (pvertex.GetRr()<0.5){
3533 // continue;
3534 //}
3535 // cout<<"vertex.GetChi2V0()"<<vertex.GetChi2V0()<<endl;
3536 if(cpa<0.9)continue;
3537 // if (vertex.GetChi2V0() > 30) continue;
3538 // cout<<"xp+xn::"<<xp<<" "<<xn<<endl;
3539 if ((xn+xp) < 0.4) continue;
3540 if (TMath::Abs(ntrk->GetD(xPrimaryVertex,yPrimaryVertex,b))<0.05)
3541 if (TMath::Abs(ptrk->GetD(xPrimaryVertex,yPrimaryVertex,b))<0.05) continue;
3542
3543 //cout<<"pass"<<endl;
3544
3545 AliKFParticle v0GammaC;
3546 v0GammaC+=(*negKF);
3547 v0GammaC+=(*posKF);
3548 v0GammaC.SetMassConstraint(0,0.001);
3549 primVtxImprovedGamma+=v0GammaC;
3550 v0GammaC.SetProductionVertex(primVtxImprovedGamma);
3551
3552
3553 curGamma=curElecNeg+curElecPos;
3554 curGammaAt=curElecNegAt+curElecPosAt;
3555
3556 // invariant mass versus pt of K0short
3557
3558 Double_t chi2V0GammaC=100000.;
3559 if( v0GammaC.GetNDF() != 0) {
3560 chi2V0GammaC = v0GammaC.GetChi2()/v0GammaC.GetNDF();
3561 }else{
3562 cout<< "ERROR::v0K0C.GetNDF()" << endl;
3563 }
3564
3565 if(chi2V0GammaC<200 &&chi2V0GammaC>0 ){
3566 if(fHistograms != NULL){
3567 fHistograms->FillHistogram("ESD_RecalculateV0_InvMass",v0GammaC.GetMass());
3568 fHistograms->FillHistogram("ESD_RecalculateV0_Pt",v0GammaC.GetPt());
3569 fHistograms->FillHistogram("ESD_RecalculateV0_E_dEdxP",curElecNegAt.P(),negTrack->GetTPCsignal());
3570 fHistograms->FillHistogram("ESD_RecalculateV0_P_dEdxP",curElecPosAt.P(),posTrack->GetTPCsignal());
3571 fHistograms->FillHistogram("ESD_RecalculateV0_cpa",cpa);
3572 fHistograms->FillHistogram("ESD_RecalculateV0_dca",dca);
9c1cb6f7 3573 fHistograms->FillHistogram("ESD_RecalculateV0_normdistP",normdistP);
3574 fHistograms->FillHistogram("ESD_RecalculateV0_normdistN",normdistN);
3575
3576 new((*fKFRecalculatedGammasTClone)[fKFRecalculatedGammasTClone->GetEntriesFast()]) AliKFParticle(v0GammaC);
3577 fElectronRecalculatedv1.push_back(iTracks);
3578 fElectronRecalculatedv2.push_back(jTracks);
3579 }
3580 }
3581 }
3582 }
3583
3584 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFRecalculatedGammasTClone->GetEntriesFast();firstGammaIndex++){
3585 for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFRecalculatedGammasTClone->GetEntriesFast();secondGammaIndex++){
3586 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFRecalculatedGammasTClone->At(firstGammaIndex);
3587 AliKFParticle * twoGammaDecayCandidateDaughter1 = (AliKFParticle *)fKFRecalculatedGammasTClone->At(secondGammaIndex);
3588
3589 if(fElectronRecalculatedv1[firstGammaIndex]==fElectronRecalculatedv1[secondGammaIndex]){
3590 continue;
3591 }
3592 if( fElectronRecalculatedv2[firstGammaIndex]==fElectronRecalculatedv2[secondGammaIndex]){
3593 continue;
3594 }
3595
3596 AliKFParticle twoGammaCandidate(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
3597 fHistograms->FillHistogram("ESD_RecalculateGG_InvMass",twoGammaCandidate.GetMass());
3598 fHistograms->FillHistogram("ESD_RecalculateGG_InvMass_vs_Pt",twoGammaCandidate.GetMass(),twoGammaCandidate.GetPt());
3599 }
3600 }
3601}
01b7fdcc 3602void AliAnalysisTaskGammaConversion::CalculateJetCone(Int_t gammaIndex){
6c84d371 3603 // CaculateJetCone
a0b94e5c 3604
d7d7e825 3605 Double_t cone;
3606 Double_t coneSize=0.3;
3607 Double_t ptJet=0;
a0b94e5c 3608
6c84d371 3609 // AliKFParticle * currentGamma = &fKFReconstructedGammas[gammaIndex];
3610 AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
98778c17 3611
01b7fdcc 3612 TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
a0b94e5c 3613
6c84d371 3614 AliESDtrack* leadingCharged = (AliESDtrack*)(fChargedParticles->At(fLeadingChargedIndex));
98778c17 3615
d7d7e825 3616 Double_t momLeadingCharged[3];
3617 leadingCharged->GetConstrainedPxPyPz(momLeadingCharged);
a0b94e5c 3618
01b7fdcc 3619 TVector3 momentumVectorLeadingCharged(momLeadingCharged[0],momLeadingCharged[1],momLeadingCharged[2]);
a0b94e5c 3620
01b7fdcc 3621 Double_t phi1=momentumVectorLeadingCharged.Phi();
3622 Double_t eta1=momentumVectorLeadingCharged.Eta();
3623 Double_t phi3=momentumVectorCurrentGamma.Phi();
a0b94e5c 3624
6c84d371 3625 for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
3626 AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
d7d7e825 3627 Int_t chId = fChargedParticlesId[iCh];
3628 if(fLeadingChargedIndex==chId || fLeadingChargedIndex==chId) continue;
3629 Double_t mom[3];
3630 curTrack->GetConstrainedPxPyPz(mom);
01b7fdcc 3631 TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
3632 Double_t phi2=momentumVectorChargedParticle.Phi();
3633 Double_t eta2=momentumVectorChargedParticle.Eta();
a0b94e5c 3634
3635
d7d7e825 3636 cone=100.;
3637 if( TMath::Abs(phi2 - phi1) <= ( TMath::TwoPi()-coneSize) ){
3638 cone = TMath::Sqrt( TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-phi1),2) );
3639 }else{
3640 if( (phi2 - phi1)> TMath::TwoPi()-coneSize ){
3641 cone = TMath::Sqrt( TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-TMath::TwoPi()-phi1),2) );
3642 }
3643 if( (phi2 - phi1)< -(TMath::TwoPi()-coneSize) ){
3644 cone = TMath::Sqrt( TMath::Power((eta2-eta1),2)+ TMath::Power((phi2+TMath::TwoPi()-phi1),2) );
3645 }
3646 }
a0b94e5c 3647
01b7fdcc 3648 if(cone <coneSize&& momentumVectorChargedParticle.Pt()>fMinPtJetCone ){
3649 ptJet+= momentumVectorChargedParticle.Pt();
3650 Double_t ffzHdrGam = momentumVectorChargedParticle.Pt()/momentumVectorCurrentGamma.Pt();
3651 Double_t imbalanceHdrGam=-momentumVectorChargedParticle.Dot(momentumVectorCurrentGamma)/momentumVectorCurrentGamma.Mag2();
d7d7e825 3652 fHistograms->FillHistogram("ESD_FFzHdrGam",ffzHdrGam);
3653 fHistograms->FillHistogram("ESD_ImbalanceHdrGam",imbalanceHdrGam);
a0b94e5c 3654
d7d7e825 3655 }
a0b94e5c 3656
d7d7e825 3657 Double_t dphiHdrGam=phi3-phi2;
3658 if ( dphiHdrGam < (-TMath::PiOver2())){
3659 dphiHdrGam+=(TMath::TwoPi());
3660 }
a0b94e5c 3661
d7d7e825 3662 if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
3663 dphiHdrGam-=(TMath::TwoPi());
3664 }
a0b94e5c 3665
01b7fdcc 3666 if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
d7d7e825 3667 fHistograms->FillHistogram("ESD_dphiHdrGamIsolated",dphiHdrGam);
3668 }
3669 }//track loop
a0b94e5c 3670
3671
d7d7e825 3672}
3673
3674Double_t AliAnalysisTaskGammaConversion::GetMinimumDistanceToCharge(Int_t indexHighestPtGamma){
3675 // GetMinimumDistanceToCharge
a0b94e5c 3676
d7d7e825 3677 Double_t fIsoMin=100.;
3678 Double_t ptLeadingCharged=-1.;
98778c17 3679
3680 fLeadingChargedIndex=-1;
a0b94e5c 3681
6c84d371 3682 AliKFParticle * gammaHighestPt = (AliKFParticle*)fKFReconstructedGammasTClone->At(indexHighestPtGamma);
01b7fdcc 3683 TVector3 momentumVectorgammaHighestPt(gammaHighestPt->GetPx(),gammaHighestPt->GetPy(),gammaHighestPt->GetPz());
a0b94e5c 3684
01b7fdcc 3685 Double_t phi1=momentumVectorgammaHighestPt.Phi();
3686 Double_t eta1=momentumVectorgammaHighestPt.Eta();
a0b94e5c 3687
6c84d371 3688 for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
3689 AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
d7d7e825 3690 Int_t chId = fChargedParticlesId[iCh];
3691 if(fElectronv1[indexHighestPtGamma]==chId || fElectronv2[indexHighestPtGamma]==chId) continue;
3692 Double_t mom[3];
3693 curTrack->GetConstrainedPxPyPz(mom);
01b7fdcc 3694 TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
3695 Double_t phi2=momentumVectorChargedParticle.Phi();
3696 Double_t eta2=momentumVectorChargedParticle.Eta();
d7d7e825 3697 Double_t iso=pow( (pow( (eta1-eta2),2)+ pow((phi1-phi2),2)),0.5 );
a0b94e5c 3698
01b7fdcc 3699 if(momentumVectorChargedParticle.Pt()>fMinPtIsoCone ){
d7d7e825 3700 if (iso<fIsoMin){
3701 fIsoMin=iso;
3702 }
3703 }
a0b94e5c 3704
d7d7e825 3705 Double_t dphiHdrGam=phi1-phi2;
3706 if ( dphiHdrGam < (-TMath::PiOver2())){
3707 dphiHdrGam+=(TMath::TwoPi());
3708 }
a0b94e5c 3709
d7d7e825 3710 if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
3711 dphiHdrGam-=(TMath::TwoPi());
3712 }
01b7fdcc 3713 if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
d7d7e825 3714 fHistograms->FillHistogram("ESD_dphiHdrGam",dphiHdrGam);
3715 }
a0b94e5c 3716
d7d7e825 3717 if (dphiHdrGam>0.9*TMath::Pi() && dphiHdrGam<1.1*TMath::Pi()){
a0b94e5c 3718 if (momentumVectorChargedParticle.Pt()> ptLeadingCharged && momentumVectorChargedParticle.Pt()>0.1*momentumVectorgammaHighestPt.Pt()){
01b7fdcc 3719 ptLeadingCharged=momentumVectorChargedParticle.Pt();
d7d7e825 3720 fLeadingChargedIndex=iCh;
3721 }
3722 }
a0b94e5c 3723
d7d7e825 3724 }//track loop
3725 fHistograms->FillHistogram("ESD_MinimumIsoDistance",fIsoMin);
3726 return fIsoMin;
a0b94e5c 3727
d7d7e825 3728}
3729
a0b94e5c 3730Int_t AliAnalysisTaskGammaConversion::GetIndexHighestPtGamma(){
3731 //GetIndexHighestPtGamma
3732
d7d7e825 3733 Int_t indexHighestPtGamma=-1;
01b7fdcc 3734 //Double_t
3735 fGammaPtHighest = -100.;
a0b94e5c 3736
6c84d371 3737 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
3738 AliKFParticle * gammaHighestPtCandidate = (AliKFParticle*)fKFReconstructedGammasTClone->At(firstGammaIndex);
a0b94e5c 3739 TVector3 momentumVectorgammaHighestPtCandidate(gammaHighestPtCandidate->GetPx(),gammaHighestPtCandidate->GetPy(),gammaHighestPtCandidate->GetPz());
3740 if (momentumVectorgammaHighestPtCandidate.Pt() > fGammaPtHighest){
3741 fGammaPtHighest=momentumVectorgammaHighestPtCandidate.Pt();
3742 //gammaHighestPt = gammaHighestPtCandidate;
3743 indexHighestPtGamma=firstGammaIndex;
3744 }
d7d7e825 3745 }
a0b94e5c 3746
d7d7e825 3747 return indexHighestPtGamma;
a0b94e5c 3748
d7d7e825 3749}
3750
3751
3752void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)
3753{
3754 // Terminate analysis
3755 //
3756 AliDebug(1,"Do nothing in Terminate");
3757}
3758
3759void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()
3760{
332f1f44 3761
63477d82 3762 if(fKFCreateAOD) {
3763
3764 //AOD
3765 if(!fAODGamma) fAODGamma = new TClonesArray(fOutputAODClassName.Data(), 0);
3766 else fAODGamma->Delete();
3767 fAODGamma->SetName(Form("%s_gamma", fAODBranchName.Data()));
3768
3769 if(!fAODPi0) fAODPi0 = new TClonesArray(fOutputAODClassName.Data(), 0);
3770 else fAODPi0->Delete();
3771 fAODPi0->SetName(Form("%s_Pi0", fAODBranchName.Data()));
3772
3773 if(!fAODOmega) fAODOmega = new TClonesArray(fOutputAODClassName.Data(), 0);
3774 else fAODOmega->Delete();
3775 fAODOmega->SetName(Form("%s_Omega", fAODBranchName.Data()));
3776
3777 //If delta AOD file name set, add in separate file. Else add in standard aod file.
3778 if(GetDeltaAODFileName().Length() > 0) {
3779 AddAODBranch("TClonesArray", &fAODGamma, GetDeltaAODFileName().Data());
3780 AddAODBranch("TClonesArray", &fAODPi0, GetDeltaAODFileName().Data());
3781 AddAODBranch("TClonesArray", &fAODOmega, GetDeltaAODFileName().Data());
3782 AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(GetDeltaAODFileName().Data());
3783 } else {
3784 AddAODBranch("TClonesArray", &fAODGamma);
3785 AddAODBranch("TClonesArray", &fAODPi0);
3786 AddAODBranch("TClonesArray", &fAODOmega);
3787 }
332f1f44 3788 }
04bf4381 3789
d7d7e825 3790 // Create the output container
3791 if(fOutputContainer != NULL){
3792 delete fOutputContainer;
3793 fOutputContainer = NULL;
3794 }
3795 if(fOutputContainer == NULL){
3796 fOutputContainer = new TList();
b58d3c74 3797 fOutputContainer->SetOwner(kTRUE);
d7d7e825 3798 }
3799
3800 //Adding the histograms to the output container
3801 fHistograms->GetOutputContainer(fOutputContainer);
3802
3803
3804 if(fWriteNtuple){
3805 if(fGammaNtuple == NULL){
3806 fGammaNtuple = new TNtuple("V0ntuple","V0ntuple","OnTheFly:HasVertex:NegPIDProb:PosPIDProb:X:Y:Z:R:MotherCandidateNDF:MotherCandidateChi2:MotherCandidateEnergy:MotherCandidateEta:MotherCandidatePt:MotherCandidateMass:MotherCandidateWidth:MCMotherCandidatePT:EPOpeningAngle:ElectronEnergy:ElectronPt:ElectronEta:ElectronPhi:PositronEnergy:PositronPt:PositronEta:PositronPhi:HasSameMCMother:MotherMCParticlePIDCode",50000);
3807 }
3808 if(fNeutralMesonNtuple == NULL){
3809 fNeutralMesonNtuple = new TNtuple("NeutralMesonNtuple","NeutralMesonNtuple","test");
3810 }
3811 TList * ntupleTList = new TList();
b58d3c74 3812 ntupleTList->SetOwner(kTRUE);
d7d7e825 3813 ntupleTList->SetName("Ntuple");
3814 ntupleTList->Add((TNtuple*)fGammaNtuple);
3815 fOutputContainer->Add(ntupleTList);
3816 }
3817
3818 fOutputContainer->SetName(GetName());
3819}
3820
77ac6f3e 3821Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(const TParticle* const daughter0, const TParticle* const daughter1) const{
d7d7e825 3822 //helper function
3823 TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());
3824 TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());
3825 return v3D0.Angle(v3D1);
3826}
3827
3828void AliAnalysisTaskGammaConversion::CheckV0Efficiency(){
3829 // see header file for documentation
5e55d806 3830
d7d7e825 3831 vector<Int_t> indexOfGammaParticle;
a0b94e5c 3832
d7d7e825 3833 fStack = fV0Reader->GetMCStack();
a0b94e5c 3834
d7d7e825 3835 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
3836 return; // aborts if the primary vertex does not have contributors.
3837 }
a0b94e5c 3838
d7d7e825 3839 for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {
3840 TParticle* particle = (TParticle *)fStack->Particle(iTracks);
3841 if(particle->GetPdgCode()==22){ //Gamma
3842 if(particle->GetNDaughters() >= 2){
3843 TParticle* electron=NULL;
3844 TParticle* positron=NULL;
3845 for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
3846 TParticle *tmpDaughter = fStack->Particle(daughterIndex);
3847 if(tmpDaughter->GetUniqueID() == 5){
3848 if(tmpDaughter->GetPdgCode() == 11){
3849 electron = tmpDaughter;
3850 }
3851 else if(tmpDaughter->GetPdgCode() == -11){
3852 positron = tmpDaughter;
3853 }
3854 }
3855 }
3856 if(electron!=NULL && positron!=0){
3857 if(electron->R()<160){
3858 indexOfGammaParticle.push_back(iTracks);
3859 }
3860 }
3861 }
3862 }
3863 }
a0b94e5c 3864
d7d7e825 3865 Int_t nFoundGammas=0;
3866 Int_t nNotFoundGammas=0;
a0b94e5c 3867
d7d7e825 3868 Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
3869 for(Int_t i=0;i<numberOfV0s;i++){
3870 fV0Reader->GetV0(i);
a0b94e5c 3871
d7d7e825 3872 if(fV0Reader->HasSameMCMother() == kFALSE){
3873 continue;
3874 }
a0b94e5c 3875
d7d7e825 3876 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
3877 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
a0b94e5c 3878
d7d7e825 3879 if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
3880 continue;
3881 }
3882 if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
3883 continue;
3884 }
a0b94e5c 3885
d7d7e825 3886 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
3887 //TParticle * v0Gamma = fV0Reader->GetMotherMCParticle();
3888 for(UInt_t mcIndex=0;mcIndex<indexOfGammaParticle.size();mcIndex++){
3889 if(negativeMC->GetFirstMother()==indexOfGammaParticle[mcIndex]){
3890 nFoundGammas++;
3891 }
3892 else{
3893 nNotFoundGammas++;
3894 }
3895 }
3896 }
3897 }
d7d7e825 3898}
3899
3900
3901void AliAnalysisTaskGammaConversion::ProcessGammaElectronsForChicAnalysis(){
3902 // see header file for documantation
a0b94e5c 3903
d7d7e825 3904 fESDEvent = fV0Reader->GetESDEvent();
a0b94e5c 3905
3906
3907 TClonesArray * vESDeNegTemp = new TClonesArray("AliESDtrack",0);
6c84d371 3908 TClonesArray * vESDePosTemp = new TClonesArray("AliESDtrack",0);
3909 TClonesArray * vESDxNegTemp = new TClonesArray("AliESDtrack",0);
3910 TClonesArray * vESDxPosTemp = new TClonesArray("AliESDtrack",0);
3911 TClonesArray * vESDeNegNoJPsi = new TClonesArray("AliESDtrack",0);
3912 TClonesArray * vESDePosNoJPsi = new TClonesArray("AliESDtrack",0);
a0b94e5c 3913
6c84d371 3914 /*
3915 vector <AliESDtrack*> vESDeNegTemp(0);
3916 vector <AliESDtrack*> vESDePosTemp(0);
3917 vector <AliESDtrack*> vESDxNegTemp(0);
3918 vector <AliESDtrack*> vESDxPosTemp(0);
3919 vector <AliESDtrack*> vESDeNegNoJPsi(0);
3920 vector <AliESDtrack*> vESDePosNoJPsi(0);
3921 */
a0b94e5c 3922
3923
d7d7e825 3924 fHistograms->FillTable("Table_Electrons",0);//Count number of Events
a0b94e5c 3925
d7d7e825 3926 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
3927 AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
a0b94e5c 3928
d7d7e825 3929 if(!curTrack){
3930 //print warning here
3931 continue;
3932 }
a0b94e5c 3933
d7d7e825 3934 double p[3];if(!curTrack->GetConstrainedPxPyPz(p))continue;
3935 double r[3];curTrack->GetConstrainedXYZ(r);
a0b94e5c 3936
d7d7e825 3937 TVector3 rXYZ(r);
a0b94e5c 3938
d7d7e825 3939 fHistograms->FillTable("Table_Electrons",4);//Count number of ESD tracks
a0b94e5c 3940
d7d7e825 3941 Bool_t flagKink = kTRUE;
3942 Bool_t flagTPCrefit = kTRUE;
3943 Bool_t flagTRDrefit = kTRUE;
3944 Bool_t flagITSrefit = kTRUE;
3945 Bool_t flagTRDout = kTRUE;
3946 Bool_t flagVertex = kTRUE;
a0b94e5c 3947
3948
d7d7e825 3949 //Cuts ---------------------------------------------------------------
a0b94e5c 3950
d7d7e825 3951 if(curTrack->GetKinkIndex(0) > 0){
3952 fHistograms->FillHistogram("Table_Electrons",5);//Count kink
3953 flagKink = kFALSE;
3954 }
a0b94e5c 3955
d7d7e825 3956 ULong_t trkStatus = curTrack->GetStatus();
a0b94e5c 3957
d7d7e825 3958 ULong_t tpcRefit = (trkStatus & AliESDtrack::kTPCrefit);
a0b94e5c 3959
d7d7e825 3960 if(!tpcRefit){
3961 fHistograms->FillHistogram("Table_Electrons",9);//Count not TPCrefit
3962 flagTPCrefit = kFALSE;
3963 }
a0b94e5c 3964
d7d7e825 3965 ULong_t itsRefit = (trkStatus & AliESDtrack::kITSrefit);
3966 if(!itsRefit){
3967 fHistograms->FillHistogram("Table_Electrons",10);//Count not ITSrefit
3968 flagITSrefit = kFALSE;
3969 }
a0b94e5c 3970
d7d7e825 3971 ULong_t trdRefit = (trkStatus & AliESDtrack::kTRDrefit);
a0b94e5c 3972
d7d7e825 3973 if(!trdRefit){
3974 fHistograms->FillHistogram("Table_Electrons",8); //Count not TRDrefit
3975 flagTRDrefit = kFALSE;
3976 }
a0b94e5c 3977
d7d7e825 3978 ULong_t trdOut = (trkStatus & AliESDtrack::kTRDout);
a0b94e5c 3979
d7d7e825 3980 if(!trdOut) {
3981 fHistograms->FillHistogram("Table_Electrons",7); //Count not TRDout
3982 flagTRDout = kFALSE;
3983 }
a0b94e5c 3984
d7d7e825 3985 double nsigmaToVxt = GetSigmaToVertex(curTrack);
a0b94e5c 3986
d7d7e825 3987 if(nsigmaToVxt > 3){
3988 fHistograms->FillHistogram("Table_Electrons",6); //Count Tracks with number of sigmas > 3
3989 flagVertex = kFALSE;
3990 }
a0b94e5c 3991
d7d7e825 3992 if(! (flagKink && flagTPCrefit && flagITSrefit && flagTRDrefit && flagTRDout && flagVertex ) ) continue;
3993 fHistograms->FillHistogram("Table_Electrons",11);//Count Tracks passed Cuts
a0b94e5c 3994
3995
d7d7e825 3996 Stat_t pid, weight;
3997 GetPID(curTrack, pid, weight);
a0b94e5c 3998
d7d7e825 3999 if(pid!=0){
4000 fHistograms->FillHistogram("Table_Electrons",12); //Count Tracks with pid != 0
4001 }
a0b94e5c 4002
d7d7e825 4003 if(pid == 0){
4004 fHistograms->FillHistogram("Table_Electrons",13); //Count Tracks with pid != 0
4005 }
a0b94e5c 4006
4007
4008
4009
a0b94e5c 4010
4011
d7d7e825 4012 TLorentzVector curElec;
4013 curElec.SetXYZM(p[0],p[1],p[2],fElectronMass);
a0b94e5c 4014
4015
113d8432 4016 if(fDoMCTruth){
4017 Int_t labelMC = TMath::Abs(curTrack->GetLabel());
4018 TParticle* curParticle = fStack->Particle(labelMC);
4019 if(curTrack->GetSign() > 0){
4020 if( pid == 0){
4021 fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
4022 fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
4023 }
4024 else{
4025 fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
4026 fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
4027 }
4028 }
4029 }
a0b94e5c 4030
4031
d7d7e825 4032 if(curTrack->GetSign() > 0){
a0b94e5c 4033
6c84d371 4034 // vESDxPosTemp.push_back(curTrack);
4035 new((*vESDxPosTemp)[vESDxPosTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
a0b94e5c 4036
d7d7e825 4037 if( pid == 0){
a0b94e5c 4038
d7d7e825 4039 fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
4040 fHistograms->FillHistogram("ESD_ElectronPosPt",curElec.Pt());
113d8432 4041 // fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
d7d7e825 4042 fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
113d8432 4043 // fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
6c84d371 4044 // vESDePosTemp.push_back(curTrack);
4045 new((*vESDePosTemp)[vESDePosTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
a0b94e5c 4046
d7d7e825 4047 }
a0b94e5c 4048
d7d7e825 4049 }
4050 else {
5e55d806 4051
6c84d371 4052 new((*vESDxNegTemp)[vESDxNegTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
a0b94e5c 4053
d7d7e825 4054 if( pid == 0){
a0b94e5c 4055
d7d7e825 4056 fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
4057 fHistograms->FillHistogram("ESD_ElectronNegPt",curElec.Pt());
d7d7e825 4058 fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
6c84d371 4059 new((*vESDeNegTemp)[vESDeNegTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
a0b94e5c 4060
d7d7e825 4061 }
a0b94e5c 4062
d7d7e825 4063 }
a0b94e5c 4064
d7d7e825 4065 }
a0b94e5c 4066
4067
d7d7e825 4068 Bool_t ePosJPsi = kFALSE;
4069 Bool_t eNegJPsi = kFALSE;
4070 Bool_t ePosPi0 = kFALSE;
4071 Bool_t eNegPi0 = kFALSE;
4072
4073 UInt_t iePosJPsi=0,ieNegJPsi=0,iePosPi0=0,ieNegPi0=0;
a0b94e5c 4074
6c84d371 4075 for(Int_t iNeg=0; iNeg < vESDeNegTemp->GetEntriesFast(); iNeg++){
4076 if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetPdgCode() == 11)
4077 if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetMother(0) > -1){
4078 Int_t labelMother = fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetMother(0);
d7d7e825 4079 TParticle* partMother = fStack ->Particle(labelMother);
4080 if (partMother->GetPdgCode() == 111){
4081 ieNegPi0 = iNeg;
4082 eNegPi0 = kTRUE;
4083 }
4084 if(partMother->GetPdgCode() == 443){ //Mother JPsi
4085 fHistograms->FillTable("Table_Electrons",14);
4086 ieNegJPsi = iNeg;
4087 eNegJPsi = kTRUE;
4088 }
4089 else{
6c84d371 4090 // vESDeNegNoJPsi.push_back(vESDeNegTemp[iNeg]);
4091 new((*vESDeNegNoJPsi)[vESDeNegNoJPsi->GetEntriesFast()]) AliESDtrack(*(AliESDtrack*)(vESDeNegTemp->At(iNeg)));
d7d7e825 4092 // cout<<"ESD No Positivo JPsi "<<endl;
4093 }
a0b94e5c 4094
d7d7e825 4095 }
4096 }
a0b94e5c 4097
6c84d371 4098 for(Int_t iPos=0; iPos < vESDePosTemp->GetEntriesFast(); iPos++){
4099 if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetPdgCode() == -11)
4100 if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetMother(0) > -1){
4101 Int_t labelMother = fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetMother(0);
d7d7e825 4102 TParticle* partMother = fStack ->Particle(labelMother);
4103 if (partMother->GetPdgCode() == 111){
4104 iePosPi0 = iPos;
4105 ePosPi0 = kTRUE;
4106 }
4107 if(partMother->GetPdgCode() == 443){ //Mother JPsi
4108 fHistograms->FillTable("Table_Electrons",15);
4109 iePosJPsi = iPos;
4110 ePosJPsi = kTRUE;
4111 }
4112 else{
6c84d371 4113 // vESDePosNoJPsi.push_back(vESDePosTemp[iPos]);
4114 new((*vESDePosNoJPsi)[vESDePosNoJPsi->GetEntriesFast()]) AliESDtrack(*(AliESDtrack*)(vESDePosTemp->At(iPos)));
d7d7e825 4115 // cout<<"ESD No Negativo JPsi "<<endl;
4116 }
a0b94e5c 4117
d7d7e825 4118 }
4119 }
4120
4121 if( eNegJPsi && ePosJPsi ){
4122 TVector3 tempeNegV,tempePosV;
6c84d371 4123 tempeNegV.SetXYZ(((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Px(),((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Py(),((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Pz());
4124 tempePosV.SetXYZ(((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Px(),((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Py(),((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Pz());
d7d7e825 4125 fHistograms->FillTable("Table_Electrons",16);
4126 fHistograms->FillHistogram("ESD_ElectronPosNegJPsiAngle",tempeNegV.Angle(tempePosV));
6c84d371 4127 fHistograms->FillHistogram("MC_ElectronPosNegJPsiAngle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->GetLabel())),
4128 fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->GetLabel()))));
d7d7e825 4129 }
4130
4131 if( eNegPi0 && ePosPi0 ){
4132 TVector3 tempeNegV,tempePosV;
6c84d371 4133 tempeNegV.SetXYZ(((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Px(),((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Py(),((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Pz());
4134 tempePosV.SetXYZ(((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Px(),((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Py(),((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Pz());
d7d7e825 4135 fHistograms->FillHistogram("ESD_ElectronPosNegPi0Angle",tempeNegV.Angle(tempePosV));
6c84d371 4136 fHistograms->FillHistogram("MC_ElectronPosNegPi0Angle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->GetLabel())),
4137 fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->GetLabel()))));
d7d7e825 4138 }
a0b94e5c 4139
4140
d7d7e825 4141 FillAngle("ESD_eNegePosAngleBeforeCut",GetTLorentzVector(vESDeNegTemp),GetTLorentzVector(vESDePosTemp));
a0b94e5c 4142
6c84d371 4143 CleanWithAngleCuts(*vESDeNegTemp,*vESDePosTemp,*fKFReconstructedGammasTClone);
d7d7e825 4144
6c84d371 4145 // vector <TLorentzVector> vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectron);
4146 // vector <TLorentzVector> vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectron);
a0b94e5c 4147
6c84d371 4148 TClonesArray vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectronTClone);
4149 TClonesArray vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectronTClone);
a0b94e5c 4150
4151
d7d7e825 4152 FillAngle("ESD_eNegePosAngleAfterCut",vCurrentTLVeNeg,vCurrentTLVePos);
a0b94e5c 4153
4154
4155
4156
d7d7e825 4157 //FillAngle("ESD_eNegePosAngleAfterCut",CurrentTLVeNeg,CurrentTLVePos);
a0b94e5c 4158
4159
d7d7e825 4160 FillElectronInvMass("ESD_InvMass_ePluseMinus",vCurrentTLVeNeg,vCurrentTLVePos);
4161 FillElectronInvMass("ESD_InvMass_xPlusxMinus",GetTLorentzVector(vESDxNegTemp),GetTLorentzVector(vESDxPosTemp));
a0b94e5c 4162
4163
4164
6c84d371 4165 FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusChiC","ESD_InvMass_GammaePluseMinusChiCDiff",*fKFReconstructedGammasCutTClone,vCurrentTLVeNeg,vCurrentTLVePos);
a0b94e5c 4166
d7d7e825 4167 FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusPi0","ESD_InvMass_GammaePluseMinusPi0Diff",
6c84d371 4168 *fKFReconstructedGammasCutTClone,vCurrentTLVeNeg,vCurrentTLVePos);
a0b94e5c 4169
d7d7e825 4170 //BackGround
a0b94e5c 4171
d7d7e825 4172 //Like Sign e+e-
4173 ElectronBackground("ESD_ENegBackground",vCurrentTLVeNeg);
4174 ElectronBackground("ESD_EPosBackground",vCurrentTLVePos);
4175 ElectronBackground("ESD_EPosENegBackground",vCurrentTLVeNeg);
4176 ElectronBackground("ESD_EPosENegBackground",vCurrentTLVePos);
a0b94e5c 4177
d7d7e825 4178 // Like Sign e+e- no JPsi
4179 ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDeNegNoJPsi));
4180 ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDePosNoJPsi));
a0b94e5c 4181
d7d7e825 4182 //Mixed Event
a0b94e5c 4183
6c84d371 4184 if( fCurrentEventPosElectronTClone->GetEntriesFast() > 0 && fCurrentEventNegElectronTClone->GetEntriesFast() > 0 && fKFReconstructedGammasCutTClone->GetEntriesFast() > 0 ){
d7d7e825 4185 FillGammaElectronInvMass("ESD_EPosENegGammaBackgroundMX","ESD_EPosENegGammaBackgroundMXDiff",
6c84d371 4186 *fKFReconstructedGammasCutTClone,*fPreviousEventTLVNegElectronTClone,*fPreviousEventTLVPosElectronTClone);
4187 *fPreviousEventTLVNegElectronTClone = vCurrentTLVeNeg;
4188 *fPreviousEventTLVPosElectronTClone = vCurrentTLVePos;
a0b94e5c 4189
d7d7e825 4190 }
a0b94e5c 4191
d7d7e825 4192 /*
4193 //Photons P
4194 Double_t vtx[3];
4195 vtx[0]=0;vtx[1]=0;vtx[2]=0;
4196 for(UInt_t i=0;i<fKFReconstructedGammasChic.size();i++){
a0b94e5c 4197
d7d7e825 4198 // if(fMCGammaChicTempCut[i]->GetMother(0) < 0) continue;
a0b94e5c 4199
4200
4201
d7d7e825 4202 Int_t tempLabel = fStack->Particle(fMCGammaChicTempCut[i]->GetMother(0))->GetPdgCode();
4203 // cout<<" Label Pedro Gonzalez " <<tempLabel <<endl;
a0b94e5c 4204
d7d7e825 4205 // cout<<" Label Distance"<<fKFReconstructedGammasChic[i].GetDistanceFromVertex(vtx)<<endl;
a0b94e5c 4206
d7d7e825 4207 if( tempLabel == 10441 || tempLabel == 20443 || tempLabel == 445 )
a0b94e5c 4208
d7d7e825 4209 fHistograms->FillHistogram("ESD_PhotonsMomentum",fKFReconstructedGammasChic[i].GetMomentum());
a0b94e5c 4210
4211
d7d7e825 4212 }
a0b94e5c 4213
4214
d7d7e825 4215 */
1e7846f4 4216
4217
4218 vESDeNegTemp->Delete();
4219 vESDePosTemp->Delete();
4220 vESDxNegTemp->Delete();
4221 vESDxPosTemp->Delete();
4222 vESDeNegNoJPsi->Delete();
4223 vESDePosNoJPsi->Delete();
4224
4225 delete vESDeNegTemp;
4226 delete vESDePosTemp;
4227 delete vESDxNegTemp;
4228 delete vESDxPosTemp;
4229 delete vESDeNegNoJPsi;
4230 delete vESDePosNoJPsi;
d7d7e825 4231}
4232
6c84d371 4233/*
4234 void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,vector <TLorentzVector> tlVeNeg, vector <TLorentzVector> tlVePos){
d7d7e825 4235 //see header file for documentation
4236 for( UInt_t iNeg=0; iNeg < tlVeNeg.size(); iNeg++){
6c84d371 4237 for (UInt_t iPos=0; iPos < tlVePos.size(); iPos++){
4238 fHistograms->FillHistogram(histoName.Data(),tlVeNeg[iNeg].Vect().Angle(tlVePos[iPos].Vect()));
4239 }
4240 }
4241 }
4242*/
4243void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,TClonesArray const tlVeNeg, TClonesArray const tlVePos){
4244 //see header file for documentation
4245 for( Int_t iNeg=0; iNeg < tlVeNeg.GetEntriesFast(); iNeg++){
4246 for (Int_t iPos=0; iPos < tlVePos.GetEntriesFast(); iPos++){
4247 fHistograms->FillHistogram(histoName.Data(),((TLorentzVector*)(tlVeNeg.At(iNeg)))->Vect().Angle(((TLorentzVector*)(tlVePos.At(iPos)))->Vect()));
d7d7e825 4248 }
4249 }
4250}
6c84d371 4251void AliAnalysisTaskGammaConversion::FillElectronInvMass(TString histoName, TClonesArray const eNeg, TClonesArray const ePos){
d7d7e825 4252 //see header file for documentation
6c84d371 4253 for( Int_t n=0; n < eNeg.GetEntriesFast(); n++){
4254 TLorentzVector en = (*(TLorentzVector*)(eNeg.At(n)));
4255 for (Int_t p=0; p < ePos.GetEntriesFast(); p++){
4256 TLorentzVector ep = (*(TLorentzVector*)(ePos.At(p)));
d7d7e825 4257 TLorentzVector np = ep + en;
4258 fHistograms->FillHistogram(histoName.Data(),np.M());
4259 }
4260 }
d7d7e825 4261}
4262
6c84d371 4263void AliAnalysisTaskGammaConversion::FillGammaElectronInvMass(TString histoMass,TString histoDiff,TClonesArray const fKFGammas,
4264 TClonesArray const tlVeNeg,TClonesArray const tlVePos)
d7d7e825 4265{
4266 //see header file for documentation
a0b94e5c 4267
6c84d371 4268 for( Int_t iNeg=0; iNeg < tlVeNeg.GetEntriesFast(); iNeg++ ){
a0b94e5c 4269
6c84d371 4270 for (Int_t iPos=0; iPos < tlVePos.GetEntriesFast(); iPos++){
a0b94e5c 4271
6c84d371 4272 TLorentzVector xy = *((TLorentzVector *)(tlVePos.At(iPos))) + *((TLorentzVector *)(tlVeNeg.At(iNeg)));
a0b94e5c 4273
6c84d371 4274 for (Int_t iGam=0; iGam < fKFGammas.GetEntriesFast(); iGam++){
a0b94e5c 4275
6c84d371 4276 // AliKFParticle * gammaCandidate = &fKFGammas[iGam];
4277 AliKFParticle * gammaCandidate = (AliKFParticle *)(fKFGammas.At(iGam));
d7d7e825 4278 TLorentzVector g;
a0b94e5c 4279
d7d7e825 4280 g.SetXYZM(gammaCandidate->GetPx(),gammaCandidate->GetPy(),gammaCandidate->GetPz(),fGammaMass);
4281 TLorentzVector xyg = xy + g;
4282 fHistograms->FillHistogram(histoMass.Data(),xyg.M());
4283 fHistograms->FillHistogram(histoDiff.Data(),(xyg.M()-xy.M()));
4284 }
4285 }
4286 }
a0b94e5c 4287
d7d7e825 4288}
6c84d371 4289void AliAnalysisTaskGammaConversion::ElectronBackground(TString hBg, TClonesArray e)
d7d7e825 4290{
4291 // see header file for documentation
6c84d371 4292 for(Int_t i=0; i < e.GetEntriesFast(); i++)
d7d7e825 4293 {
6c84d371 4294 for (Int_t j=i+1; j < e.GetEntriesFast(); j++)
d7d7e825 4295 {
6c84d371 4296 TLorentzVector ee = (*(TLorentzVector*)(e.At(i))) + (*(TLorentzVector*)(e.At(j)));
a0b94e5c 4297
d7d7e825 4298 fHistograms->FillHistogram(hBg.Data(),ee.M());
4299 }
4300 }
4301}
4302
4303
6c84d371 4304void AliAnalysisTaskGammaConversion::CleanWithAngleCuts(TClonesArray const negativeElectrons,
4305 TClonesArray const positiveElectrons,
4306 TClonesArray const gammas){
d7d7e825 4307 // see header file for documentation
a0b94e5c 4308
6c84d371 4309 UInt_t sizeN = negativeElectrons.GetEntriesFast();
4310 UInt_t sizeP = positiveElectrons.GetEntriesFast();
4311 UInt_t sizeG = gammas.GetEntriesFast();
a0b94e5c 4312
4313
4314
d7d7e825 4315 vector <Bool_t> xNegBand(sizeN);
4316 vector <Bool_t> xPosBand(sizeP);
4317 vector <Bool_t> gammaBand(sizeG);
a0b94e5c 4318
4319
d7d7e825 4320 for(UInt_t iNeg=0; iNeg < sizeN; iNeg++) xNegBand[iNeg]=kTRUE;
4321 for(UInt_t iPos=0; iPos < sizeP; iPos++) xPosBand[iPos]=kTRUE;
4322 for(UInt_t iGam=0; iGam < sizeG; iGam++) gammaBand[iGam]=kTRUE;
d7d7e825 4323
a0b94e5c 4324
4325 for(UInt_t iPos=0; iPos < sizeP; iPos++){
4326
6c84d371 4327 Double_t aP[3];
4328 ((AliESDtrack*)(positiveElectrons.At(iPos)))->GetConstrainedPxPyPz(aP);
a0b94e5c 4329
d7d7e825 4330 TVector3 ePosV(aP[0],aP[1],aP[2]);
a0b94e5c 4331
d7d7e825 4332 for(UInt_t iNeg=0; iNeg < sizeN; iNeg++){
a0b94e5c 4333
6c84d371 4334 Double_t aN[3];
4335 ((AliESDtrack*)(negativeElectrons.At(iNeg)))->GetConstrainedPxPyPz(aN);
d7d7e825 4336 TVector3 eNegV(aN[0],aN[1],aN[2]);
a0b94e5c 4337
d7d7e825 4338 if(ePosV.Angle(eNegV) < 0.05){ //e+e- from gamma
4339 xPosBand[iPos]=kFALSE;
4340 xNegBand[iNeg]=kFALSE;
4341 }
a0b94e5c 4342
d7d7e825 4343 for(UInt_t iGam=0; iGam < sizeG; iGam++){
6c84d371 4344 AliKFParticle* gammaCandidate = (AliKFParticle*)gammas.At(iGam);
d7d7e825 4345 TVector3 gammaCandidateVector(gammaCandidate->Px(),gammaCandidate->Py(),gammaCandidate->Pz());
4346 if(ePosV.Angle(gammaCandidateVector) < 0.05 || eNegV.Angle(gammaCandidateVector) < 0.05)
4347 gammaBand[iGam]=kFALSE;
4348 }
4349 }
4350 }
a0b94e5c 4351
4352
4353
4354
d7d7e825 4355 for(UInt_t iPos=0; iPos < sizeP; iPos++){
4356 if(xPosBand[iPos]){
6c84d371 4357 new((*fCurrentEventPosElectronTClone)[fCurrentEventPosElectronTClone->GetEntriesFast()]) AliESDtrack((*(AliESDtrack*)(positiveElectrons.At(iPos))));
4358 // fCurrentEventPosElectron.push_back(positiveElectrons[iPos]);
d7d7e825 4359 }
4360 }
4361 for(UInt_t iNeg=0;iNeg < sizeN; iNeg++){
4362 if(xNegBand[iNeg]){
6c84d371 4363 new((*fCurrentEventNegElectronTClone)[fCurrentEventNegElectronTClone->GetEntriesFast()]) AliESDtrack((*(AliESDtrack*)(negativeElectrons.At(iNeg))));
4364 // fCurrentEventNegElectron.push_back(negativeElectrons[iNeg]);
d7d7e825 4365 }
4366 }
4367 for(UInt_t iGam=0; iGam < sizeG; iGam++){
4368 if(gammaBand[iGam]){
6c84d371 4369 new((*fKFReconstructedGammasCutTClone)[fKFReconstructedGammasCutTClone->GetEntriesFast()]) AliKFParticle((*(AliKFParticle*)(gammas.At(iGam))));
4370 //fKFReconstructedGammasCut.push_back(*(AliKFParticle*)gammas->At(iGam));
d7d7e825 4371 }
4372 }
4373}
4374
4375
77ac6f3e 4376void AliAnalysisTaskGammaConversion::GetPID(const AliESDtrack *track, Stat_t &pid, Stat_t &weight)
d7d7e825 4377{
4378 // see header file for documentation
4379 pid = -1;
4380 weight = -1;
a0b94e5c 4381
d7d7e825 4382 double wpart[5];
4383 double wpartbayes[5];
a0b94e5c 4384
d7d7e825 4385 //get probability of the diffenrent particle types
4386 track->GetESDpid(wpart);
a0b94e5c 4387
d7d7e825 4388 // Tentative particle type "concentrations"
4389 double c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
a0b94e5c 4390
d7d7e825 4391 //Bayes' formula
4392 double rcc = 0.;
4393 for (int i = 0; i < 5; i++)
4394 {
4395 rcc+=(c[i] * wpart[i]);
4396 }
a0b94e5c 4397
4398
4399
d7d7e825 4400 for (int i=0; i<5; i++) {
4a6157dc 4401 if( rcc>0 || rcc<0){//Kenneth: not sure if the rcc<0 should be there, this is from fixing a coding violation where we are not allowed to say: rcc!=0 (RC19)
d7d7e825 4402 wpartbayes[i] = c[i] * wpart[i] / rcc;
4403 }
4404 }
a0b94e5c 4405
4406
4407
d7d7e825 4408 Float_t max=0.;
4409 int ipid=-1;
4410 //find most probable particle in ESD pid
4411 //0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton
4412 for (int i = 0; i < 5; i++)
4413 {
4414 if (wpartbayes[i] > max)
4415 {
a0b94e5c 4416 ipid = i;
4417 max = wpartbayes[i];
d7d7e825 4418 }
4419 }
a0b94e5c 4420
d7d7e825 4421 pid = ipid;
4422 weight = max;
4423}
77ac6f3e 4424double AliAnalysisTaskGammaConversion::GetSigmaToVertex(const AliESDtrack* t)
d7d7e825 4425{
4426 // Calculates the number of sigma to the vertex.
a0b94e5c 4427
d7d7e825 4428 Float_t b[2];
4429 Float_t bRes[2];
4430 Float_t bCov[3];
4431 t->GetImpactParameters(b,bCov);
4432 if (bCov[0]<=0 || bCov[2]<=0) {
4433 AliDebug(1, "Estimated b resolution lower or equal zero!");
4434 bCov[0]=0; bCov[2]=0;
4435 }
4436 bRes[0] = TMath::Sqrt(bCov[0]);
4437 bRes[1] = TMath::Sqrt(bCov[2]);
a0b94e5c 4438
d7d7e825 4439 // -----------------------------------
4440 // How to get to a n-sigma cut?
4441 //
4442 // The accumulated statistics from 0 to d is
4443 //
4444 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
4445 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
4446 //
4447 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
4448 // Can this be expressed in a different way?
a0b94e5c 4449
d7d7e825 4450 if (bRes[0] == 0 || bRes[1] ==0)
4451 return -1;
a0b94e5c 4452
d7d7e825 4453 double d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
a0b94e5c 4454
d7d7e825 4455 // stupid rounding problem screws up everything:
4456 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
4457 if (TMath::Exp(-d * d / 2) < 1e-10)
4458 return 1000;
a0b94e5c 4459
4460
d7d7e825 4461 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
4462 return d;
4463}
6c84d371 4464
4465//vector <TLorentzVector> AliAnalysisTaskGammaConversion::GetTLorentzVector(vector <AliESDtrack*> esdTrack){
4466TClonesArray AliAnalysisTaskGammaConversion::GetTLorentzVector(TClonesArray *const esdTrack){
01b7fdcc 4467 //Return TLoresntz vector of track?
6c84d371 4468 // vector <TLorentzVector> tlVtrack(0);
4469 TClonesArray array("TLorentzVector",0);
a0b94e5c 4470
6c84d371 4471 for(Int_t itrack=0; itrack < esdTrack->GetEntriesFast(); itrack++){
4472 double p[3];
4473 //esdTrack[itrack]->GetConstrainedPxPyPz(p);
4474 ((AliESDtrack*)(esdTrack->At(itrack)))->GetConstrainedPxPyPz(p);
d7d7e825 4475 TLorentzVector currentTrack;
01b7fdcc 4476 currentTrack.SetXYZM(p[0],p[1],p[2],fElectronMass);
6c84d371 4477 new((array)[array.GetEntriesFast()]) TLorentzVector(currentTrack);
4478 // tlVtrack.push_back(currentTrack);
d7d7e825 4479 }
a0b94e5c 4480
6c84d371 4481 return array;
a0b94e5c 4482
6c84d371 4483 // return tlVtrack;
d7d7e825 4484}
77ac6f3e 4485Int_t AliAnalysisTaskGammaConversion::GetProcessType(const AliMCEvent * mcEvt) {
3c45d101 4486
4487 // Determine if the event was generated with pythia or phojet and return the process type
4488
4489 // Check if mcEvt is fine
4490 if (!mcEvt) {
4491 AliFatal("NULL mc event");
4492 }
4493
4494 // Determine if it was a pythia or phojet header, and return the correct process type
4495 AliGenPythiaEventHeader * headPy = 0;
4496 AliGenDPMjetEventHeader * headPho = 0;
4497 AliGenEventHeader * htmp = mcEvt->GenEventHeader();
4498 if(!htmp) {
4499 AliFatal("Cannot Get MC Header!!");
4500 }
4501 if( TString(htmp->IsA()->GetName()) == "AliGenPythiaEventHeader") {
4502 headPy = (AliGenPythiaEventHeader*) htmp;
4503 } else if (TString(htmp->IsA()->GetName()) == "AliGenDPMjetEventHeader") {
4504 headPho = (AliGenDPMjetEventHeader*) htmp;
4505 } else {
4506 AliError("Unknown header");
4507 }
4508
4509 // Determine process type
4510 if(headPy) {
4511 if(headPy->ProcessType() == 92 || headPy->ProcessType() == 93) {
4512 // single difractive
4513 return kProcSD;
4514 } else if (headPy->ProcessType() == 94) {
4515 // double diffractive
4516 return kProcDD;
4517 }
4518 else if(headPy->ProcessType() != 92 && headPy->ProcessType() != 93 && headPy->ProcessType() != 94) {
4519 // non difractive
4520 return kProcND;
4521 }
4522 } else if (headPho) {
4523 if(headPho->ProcessType() == 5 || headPho->ProcessType() == 6 ) {
4524 // single difractive
4525 return kProcSD;
4526 } else if (headPho->ProcessType() == 7) {
4527 // double diffractive
4528 return kProcDD;
4529 } else if(headPho->ProcessType() != 5 && headPho->ProcessType() != 6 && headPho->ProcessType() != 7 ) {
4530 // non difractive
4531 return kProcND;
4532 }
4533 }
4534
4535
4536 // no process type found?
4537 AliError(Form("Unknown header: %s", htmp->IsA()->GetName()));
4538 return kProcUnknown;
4539}
6746e1e1 4540
4541
4542Int_t AliAnalysisTaskGammaConversion::CalculateMultiplicityBin(){
4543 // Get Centrality bin
4544
4545 Int_t multiplicity = 0;
4546
4547 if ( fUseMultiplicity == 1 ) {
4548
4549 if (fMultiplicity>= 0 && fMultiplicity<= 5) multiplicity=1;
4550 if (fMultiplicity>= 6 && fMultiplicity<= 9) multiplicity=2;
4551 if (fMultiplicity>=10 && fMultiplicity<=14) multiplicity=3;
4552 if (fMultiplicity>=15 && fMultiplicity<=22) multiplicity=4;
4553 if (fMultiplicity>=23 ) multiplicity=5;
4554
4555 }
4556 return multiplicity;
4557}