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