Changes for report #74527: request to commit code under TRIGGER to handle TOF trigger...
[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;
406a77e0 332 fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
d707e3cf 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/////////////////////////////
91b98865 672 if(fDoChic) {
673 if(particle->GetPdgCode() == 443){//Is JPsi
674 if(particle->GetNDaughters()==2){
675 if(TMath::Abs(fStack->Particle(particle->GetFirstDaughter())->GetPdgCode()) == 11 &&
676 TMath::Abs(fStack->Particle(particle->GetLastDaughter())->GetPdgCode()) == 11){
677
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
91b98865 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
91b98865 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 }
d7d7e825 696 }
697 }
698 }
91b98865 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 }
d7d7e825 710 }
711 }
91b98865 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
91b98865 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
91b98865 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 }
d7d7e825 734 }
a0b94e5c 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);
9da64c46 2056 if(massTwoGammaCandidate>0.1 && massTwoGammaCandidate<0.15){
2057 fHistograms->FillHistogram("ESD_Mother_alfa_Pi0", alfa);
2058 }
6de3471d 2059 fHistograms->FillHistogram("ESD_Mother_R", spaceVectorTwoGammaCandidate.Pt()); // Pt in Space == R!!!
2060 fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), spaceVectorTwoGammaCandidate.Pt());
2061 fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());
2062 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2063 fHistograms->FillHistogram("ESD_Mother_InvMass",massTwoGammaCandidate);
9c1cb6f7 2064 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2065 }
3c45d101 2066 if(alfa<0.1){
2067 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_E_alpha",massTwoGammaCandidate ,twoGammaCandidate->GetE());
2068 }
6de3471d 2069
9640a3d1 2070
9c1cb6f7 2071 if(fCalculateBackground){
5ce758b0 2072 /* Kenneth, just for testing*/
2073 AliGammaConversionBGHandler * bgHandlerTest = fV0Reader->GetBGHandler();
2074
2075 Int_t zbin= bgHandlerTest->GetZBinIndex(fV0Reader->GetVertexZ());
2076 Int_t mbin=0;
2077 Int_t multKAA=0;
2078 if(fUseTrackMultiplicityForBG == kTRUE){
2079 multKAA=fV0Reader->CountESDTracks();
2080 mbin = bgHandlerTest->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2081 }
2082 else{// means we use #v0s for multiplicity
2083 multKAA=fV0Reader->GetNGoodV0s();
2084 mbin = bgHandlerTest->GetMultiplicityBinIndex(fV0Reader->GetNGoodV0s());
2085 }
2086 // cout<<"Filling bin number "<<zbin<<" and "<<mbin<<endl;
2087 // cout<<"Corresponding to z = "<<fV0Reader->GetVertexZ()<<" and m = "<<multKAA<<endl;
6de3471d 2088 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2089 fHistograms->FillHistogram(Form("%d%dESD_Mother_InvMass",zbin,mbin),massTwoGammaCandidate);
2090 fHistograms->FillHistogram(Form("%d%dESD_Mother_InvMass_vs_Pt",zbin,mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2091 /* end Kenneth, just for testing*/
2092 fHistograms->FillHistogram(Form("%dESD_Mother_InvMass_vs_Pt",mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2093 }
5ce758b0 2094 }
2095 /* if(fCalculateBackground){
9c1cb6f7 2096 AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
2097 Int_t mbin= bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2098 fHistograms->FillHistogram(Form("%dESD_Mother_InvMass_vs_Pt",mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
5ce758b0 2099 }*/
ebcfaa7e 2100 // if(fDoNeutralMesonV0MCCheck){
2101 if(fDoMCTruth){
037dc2db 2102 //Kenneth: Checking the eta of the gamma to check the difference between 0.9 and 1.2
2103 Int_t indexKF1 = fKFReconstructedGammasV0Index.at(firstGammaIndex);
2104 if(indexKF1<fV0Reader->GetNumberOfV0s()){
2105 fV0Reader->GetV0(indexKF1);//updates to the correct v0
2106 Double_t eta1 = fV0Reader->GetMotherCandidateEta();
2107 Bool_t isRealPi0=kFALSE;
10e3319b 2108 Bool_t isRealEta=kFALSE;
037dc2db 2109 Int_t gamma1MotherLabel=-1;
2110 if(fV0Reader->HasSameMCMother() == kTRUE){
2111 //cout<<"This v0 is a real v0!!!!"<<endl;
2112 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
2113 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
2114 if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
2115 if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){
2116 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
2117 gamma1MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
2118 }
2119 }
2120 }
2121 }
2122 Int_t indexKF2 = fKFReconstructedGammasV0Index.at(secondGammaIndex);
2123 if(indexKF1 == indexKF2){
2124 cout<<"index of the two KF particles are the same.... should not happen"<<endl;
2125 }
2126 if(indexKF2<fV0Reader->GetNumberOfV0s()){
2127 fV0Reader->GetV0(indexKF2);
2128 Double_t eta2 = fV0Reader->GetMotherCandidateEta();
2129 Int_t gamma2MotherLabel=-1;
2130 if(fV0Reader->HasSameMCMother() == kTRUE){
2131 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
2132 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
2133 if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
2134 if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){
2135 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
2136 gamma2MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
2137 }
2138 }
2139 }
2140 }
2141 if(gamma1MotherLabel>=0 && gamma1MotherLabel==gamma2MotherLabel){
2142 if(fV0Reader->CheckIfPi0IsMother(gamma1MotherLabel)){
2143 isRealPi0=kTRUE;
2144 }
10e3319b 2145 if(fV0Reader->CheckIfEtaIsMother(gamma1MotherLabel)){
2146 isRealEta=kTRUE;
2147 }
2148
037dc2db 2149 }
6de3471d 2150 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2151 if(TMath::Abs(eta1)>0.9 && TMath::Abs(eta2)>0.9){
2152 // fHistograms->FillHistogram("ESD_Mother_InvMass_1212",massTwoGammaCandidate);
2153 // fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2154 if(isRealPi0 || isRealEta){
2155 fHistograms->FillHistogram("ESD_TruePi0_InvMass_1212",massTwoGammaCandidate);
2156 fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_1212",openingAngleTwoGammaCandidate);
2157 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2158 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2159 fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
d707e3cf 2160 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2161 }
6de3471d 2162
2163 if(!isRealPi0 && !isRealEta){
2164 if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
2165 fHistograms->FillHistogram("ESD_TrueBckGG_InvMass",massTwoGammaCandidate);
2166 }else{
2167 fHistograms->FillHistogram("ESD_TrueBckCont_InvMass",massTwoGammaCandidate);
2168 }
bd6d9fa3 2169 }
bd6d9fa3 2170
6de3471d 2171 }
2172 else if(TMath::Abs(eta1)>0.9 || TMath::Abs(eta2)>0.9){
ebcfaa7e 2173 // fHistograms->FillHistogram("ESD_Mother_InvMass_0912",massTwoGammaCandidate);
2174 // fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
6de3471d 2175
2176 if(isRealPi0 || isRealEta){
2177 fHistograms->FillHistogram("ESD_TruePi0_InvMass_0912",massTwoGammaCandidate);
2178 fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_0912",openingAngleTwoGammaCandidate);
2179 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2180 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2181 fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
d707e3cf 2182 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2183 }
6de3471d 2184 if(!isRealPi0 && !isRealEta){
2185 if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
2186 fHistograms->FillHistogram("ESD_TrueBckGG_InvMass",massTwoGammaCandidate);
2187 }else{
2188 fHistograms->FillHistogram("ESD_TrueBckCont_InvMass",massTwoGammaCandidate);
2189 }
bd6d9fa3 2190 }
2191 }
6de3471d 2192 else{
ebcfaa7e 2193 // fHistograms->FillHistogram("ESD_Mother_InvMass_0909",massTwoGammaCandidate);
2194 // fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
6de3471d 2195 if(isRealPi0 || isRealEta){
2196 fHistograms->FillHistogram("ESD_TruePi0_InvMass_0909",massTwoGammaCandidate);
2197 fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_0909",openingAngleTwoGammaCandidate);
2198 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2199 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2200 fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
d707e3cf 2201 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
6de3471d 2202 if(gamma1MotherLabel > fV0Reader->GetMCStack()->GetNprimary()){
2203 fHistograms->FillHistogram("ESD_TruePi0Sec_InvMass",massTwoGammaCandidate);
2204 }
d707e3cf 2205 }
6de3471d 2206 if(!isRealPi0 && !isRealEta){
2207 if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
2208 fHistograms->FillHistogram("ESD_TrueBckGG_InvMass",massTwoGammaCandidate);
2209 }else{
2210 fHistograms->FillHistogram("ESD_TrueBckCont_InvMass",massTwoGammaCandidate);
2211 }
bd6d9fa3 2212 }
037dc2db 2213 }
2214 }
2215 }
2216 }
2217 }
6de3471d 2218 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2219 if ( TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())<0.9 && TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())<0.9 ){
2220 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_Fiducial",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2221 fHistograms->FillHistogram("ESD_Mother_InvMass_Fiducial",massTwoGammaCandidate);
2222 }
2223
2224 if(TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())>0.9 && TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())>0.9){
2225 fHistograms->FillHistogram("ESD_Mother_InvMass_1212",massTwoGammaCandidate);
2226 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2227 }
2228 else if(TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())>0.9 || TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())>0.9){
2229 fHistograms->FillHistogram("ESD_Mother_InvMass_0912",massTwoGammaCandidate);
2230 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2231 }
2232 else{
2233 fHistograms->FillHistogram("ESD_Mother_InvMass_0909",massTwoGammaCandidate);
2234 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2235 }
ebcfaa7e 2236
6de3471d 2237 Double_t lowMassPi0=0.1;
2238 Double_t highMassPi0=0.15;
2239 if (massTwoGammaCandidate > lowMassPi0 && massTwoGammaCandidate < highMassPi0 ){
2240 new((*fKFReconstructedPi0sTClone)[fKFReconstructedPi0sTClone->GetEntriesFast()]) AliKFParticle(*twoGammaCandidate);
2241 fGammav1.push_back(firstGammaIndex);
2242 fGammav2.push_back(secondGammaIndex);
2243 AddPionToAOD(twoGammaCandidate, massTwoGammaCandidate, firstGammaIndex, secondGammaIndex);
2244 }
ebcfaa7e 2245 }
6272370b 2246
2e2da371 2247 }
d707e3cf 2248 //}
6de3471d 2249 delete twoGammaCandidate;
d7d7e825 2250 }
2251 }
2252}
2253
04bf4381 2254void AliAnalysisTaskGammaConversion::AddPionToAOD(AliKFParticle * pionkf, Double_t mass, Int_t daughter1, Int_t daughter2) {
2255 //See header file for documentation
2256 AliGammaConversionAODObject pion;
2257 pion.SetPx(pionkf->Px());
2258 pion.SetPy(pionkf->Py());
2259 pion.SetPz(pionkf->Pz());
2260 pion.SetChi2(pionkf->GetChi2());
2261 pion.SetE(pionkf->E());
2262 pion.SetIMass(mass);
2263 pion.SetLabel1(daughter1);
2264 //dynamic_cast<AliGammaConversionAODObject*>(fAODBranch->At(daughter1))->SetTagged(kTRUE);
2265 pion.SetLabel2(daughter2);
2266 new((*fAODPi0)[fAODPi0->GetEntriesFast()]) AliGammaConversionAODObject(pion);
2267
2268}
2269
2270
2271
6272370b 2272void AliAnalysisTaskGammaConversion::ProcessConvPHOSGammasForNeutralMesonAnalysis(){
48682642 2273/*
6272370b 2274 // see header file for documentation
2275 // Analyse Pi0 with one photon from Phos and 1 photon from conversions
2276
2277
2278
2279 Double_t vtx[3];
2280 vtx[0] = fV0Reader->GetPrimaryVertex()->GetX();
2281 vtx[1] = fV0Reader->GetPrimaryVertex()->GetY();
2282 vtx[2] = fV0Reader->GetPrimaryVertex()->GetZ();
2283
2284
2285 // Loop over all CaloClusters and consider only the PHOS ones:
2286 AliESDCaloCluster *clu;
2287 TLorentzVector pPHOS;
2288 TLorentzVector gammaPHOS;
2289 TLorentzVector gammaGammaConv;
2290 TLorentzVector pi0GammaConvPHOS;
2291 TLorentzVector gammaGammaConvBck;
2292 TLorentzVector pi0GammaConvPHOSBck;
2293
2294
2295 for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
2296 clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
2297 if ( !clu->IsPHOS() || clu->E()<0.1 ) continue;
2298 clu ->GetMomentum(pPHOS ,vtx);
2299 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2300 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2301 gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz(),0.);
2302 gammaPHOS.SetXYZM(pPHOS.Px(),pPHOS.Py(),pPHOS.Pz(),0.);
2303 pi0GammaConvPHOS=gammaGammaConv+gammaPHOS;
2304 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS",pi0GammaConvPHOS.M());
2305 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvPHOS",pi0GammaConvPHOS.M(),pi0GammaConvPHOS.Pt());
2306
2307 TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz());
2308 TVector3 v3D1(gammaPHOS.Px(),gammaPHOS.Py(),gammaPHOS.Pz());
2309 Double_t opanConvPHOS= v3D0.Angle(v3D1);
2310 if ( opanConvPHOS < 0.35){
2311 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanLow",pi0GammaConvPHOS.M());
2312 }else{
2313 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanHigh",pi0GammaConvPHOS.M());
2314 }
2315
2316 }
2317
2318 // Now the LorentVector pPHOS is obtained and can be paired with the converted proton
2319 }
2320 //==== End of the PHOS cluster selection ============
2321 TLorentzVector pEMCAL;
2322 TLorentzVector gammaEMCAL;
2323 TLorentzVector pi0GammaConvEMCAL;
2324 TLorentzVector pi0GammaConvEMCALBck;
2325
2326 for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
2327 clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
2328 if ( !clu->IsEMCAL() || clu->E()<0.1 ) continue;
2329 if (clu->GetNCells() <= 1) continue;
2330 if ( clu->GetTOF()*1e9 < 550 || clu->GetTOF()*1e9 > 750) continue;
2331
2332 clu ->GetMomentum(pEMCAL ,vtx);
2333 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2334 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2335 gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),
2336 twoGammaDecayCandidateDaughter0->Py(),
2337 twoGammaDecayCandidateDaughter0->Pz(),0.);
2338 gammaEMCAL.SetXYZM(pEMCAL.Px(),pEMCAL.Py(),pEMCAL.Pz(),0.);
2339 pi0GammaConvEMCAL=gammaGammaConv+gammaEMCAL;
2340 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL",pi0GammaConvEMCAL.M());
2341 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL",pi0GammaConvEMCAL.M(),pi0GammaConvEMCAL.Pt());
2342 TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),
2343 twoGammaDecayCandidateDaughter0->Py(),
2344 twoGammaDecayCandidateDaughter0->Pz());
2345 TVector3 v3D1(gammaEMCAL.Px(),gammaEMCAL.Py(),gammaEMCAL.Pz());
2346
2347
2348 Double_t opanConvEMCAL= v3D0.Angle(v3D1);
2349 if ( opanConvEMCAL < 0.35){
2350 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanLow",pi0GammaConvEMCAL.M());
2351 }else{
2352 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanHigh",pi0GammaConvEMCAL.M());
2353 }
2354
2355 }
48682642 2356 if(fCalculateBackground){
2357 for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2358 AliGammaConversionKFVector * previousEventV0s = fV0Reader->GetBGGoodV0s(nEventsInBG);
2359 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2360 AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2361 gammaGammaConvBck.SetXYZM(previousGoodV0.Px(),
2362 previousGoodV0.Py(),
2363 previousGoodV0.Pz(),0.);
2364 pi0GammaConvEMCALBck=gammaGammaConvBck+gammaEMCAL;
2365 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M());
2366 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M(),
2367 pi0GammaConvEMCALBck.Pt());
2368 }
6272370b 2369 }
48682642 2370
2371 // Now the LorentVector pEMCAL is obtained and can be paired with the converted proton
2372 } // end of checking if background photons are available
2373 }
6272370b 2374 //==== End of the PHOS cluster selection ============
48682642 2375*/
6272370b 2376}
2377
5ce758b0 2378void AliAnalysisTaskGammaConversion::MoveParticleAccordingToVertex(AliKFParticle *particle,AliGammaConversionBGHandler::GammaConversionVertex *vertex){
2379 //see header file for documentation
2380
2381 Double_t dx = vertex->fX - fESDEvent->GetPrimaryVertex()->GetX();
2382 Double_t dy = vertex->fY - fESDEvent->GetPrimaryVertex()->GetY();
2383 Double_t dz = vertex->fZ - fESDEvent->GetPrimaryVertex()->GetZ();
2384
2385 // cout<<"dx, dy, dz: ["<<dx<<","<<dy<<","<<dz<<"]"<<endl;
2386 particle->X() = particle->GetX() - dx;
2387 particle->Y() = particle->GetY() - dy;
2388 particle->Z() = particle->GetZ() - dz;
2389}
2390
111d75df 2391void AliAnalysisTaskGammaConversion::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle){
2392
2393 Double_t c = cos(angle);
2394 Double_t s = sin(angle);
2395
2396 Double_t A[7][ 7];
2397 for( Int_t i=0; i<7; i++ ){
2398 for( Int_t j=0; j<7; j++){
2399 A[i][j] = 0;
2400 }
2401 }
2402 for( int i=0; i<7; i++ ){
2403 A[i][i] = 1;
2404 }
2405 A[0][0] = c; A[0][1] = s;
2406 A[1][0] = -s; A[1][1] = c;
2407 A[3][3] = c; A[3][4] = s;
2408 A[4][3] = -s; A[4][4] = c;
2409
2410 Double_t AC[7][7];
2411 Double_t Ap[7];
2412
2413 for( Int_t i=0; i<7; i++ ){
2414 Ap[i] = 0;
2415 for( Int_t k=0; k<7; k++){
2416 Ap[i]+=A[i][k] * kfParticle->GetParameter(k);
2417 }
2418 }
2419
2420 for( Int_t i=0; i<7; i++){
2421 kfParticle->Parameter(i) = Ap[i];
2422 }
2423
2424 for( Int_t i=0; i<7; i++ ){
2425 for( Int_t j=0; j<7; j++ ){
2426 AC[i][j] = 0;
2427 for( Int_t k=0; k<7; k++ ){
2428 AC[i][j]+= A[i][k] * kfParticle->GetCovariance(k,j);
2429 }
2430 }
2431 }
2432
2433 for( Int_t i=0; i<7; i++ ){
2434 for( Int_t j=0; j<=i; j++ ){
2435 Double_t xx = 0;
2436 for( Int_t k=0; k<7; k++){
2437 xx+= AC[i][k]*A[j][k];
2438 }
2439 kfParticle->Covariance(i,j) = xx;
2440 }
2441 }
2442}
2443
2444
d7d7e825 2445void AliAnalysisTaskGammaConversion::CalculateBackground(){
2446 // see header file for documentation
5e55d806 2447
2448
2449 TClonesArray * currentEventV0s = fV0Reader->GetCurrentEventGoodV0s();
2450
5ce758b0 2451 AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
2452
2453 Int_t zbin= bgHandler->GetZBinIndex(fV0Reader->GetVertexZ());
2454 Int_t mbin = 0;
2455 if(fUseTrackMultiplicityForBG == kTRUE){
2456 mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2457 }
2458 else{
2459 mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->GetNGoodV0s());
2460 }
2461
111d75df 2462 if(fDoRotation == kTRUE){
3c45d101 2463 TRandom3 *random = new TRandom3();
111d75df 2464 for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2465 AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
2466 for(Int_t iCurrent2=iCurrent+1;iCurrent2<currentEventV0s->GetEntriesFast();iCurrent2++){
2467 for(Int_t nRandom=0;nRandom<nRandomEventsForBG;nRandom++){
2468
2469 AliKFParticle currentEventGoodV02 = *(AliKFParticle *)(currentEventV0s->At(iCurrent2));
3bfbe89a 2470
2471 if(fCheckBGProbability == kTRUE){
2472 Double_t massBGprob =0.;
2473 Double_t widthBGprob = 0.;
2474 AliKFParticle *backgroundCandidateProb = new AliKFParticle(currentEventGoodV0,currentEventGoodV02);
2475 backgroundCandidateProb->GetMass(massBGprob,widthBGprob);
2476 if(massBGprob>0.1 && massBGprob<0.14){
2477 if(random->Rndm()>bgHandler->GetBGProb(zbin,mbin)){
2478 delete backgroundCandidateProb;
2479 continue;
2480 }
2481 }
2482 delete backgroundCandidateProb;
2483 }
111d75df 2484
2485 Double_t nRadiansPM = nDegreesPMBackground*TMath::Pi()/180;
2486
2487 Double_t rotationValue = random->Rndm()*2*nRadiansPM + TMath::Pi()-nRadiansPM;
5ce758b0 2488
111d75df 2489 RotateKFParticle(&currentEventGoodV02,rotationValue);
5ce758b0 2490
111d75df 2491 AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,currentEventGoodV02);
5e55d806 2492
5ce758b0 2493 Double_t massBG =0.;
2494 Double_t widthBG = 0.;
2495 Double_t chi2BG =10000.;
2496 backgroundCandidate->GetMass(massBG,widthBG);
111d75df 2497
111d75df 2498 // if(backgroundCandidate->GetNDF()>0){
2499 chi2BG = backgroundCandidate->GetChi2();
2500 if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
2501
2502 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2503 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2504
2505 Double_t openingAngleBG = currentEventGoodV0.GetAngle(currentEventGoodV02);
2506
2507 Double_t rapidity;
2508 if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
2509 else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2510
dc2883e4 2511 if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ){
2512 delete backgroundCandidate;
2513 continue; // rapidity cut
2514 }
2515
111d75df 2516
2517 Double_t alfa=0.0;
2518 if( (currentEventGoodV0.GetE()+currentEventGoodV02.GetE()) != 0){
2519 alfa=TMath::Abs((currentEventGoodV0.GetE()-currentEventGoodV02.GetE())
2520 /(currentEventGoodV0.GetE()+currentEventGoodV02.GetE()));
2521 }
2522
2523
2524 if(openingAngleBG < fMinOpeningAngleGhostCut ){
2525 delete backgroundCandidate;
2526 continue; // minimum opening angle to avoid using ghosttracks
2527 }
2528
2529 // original
111d75df 2530 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
6de3471d 2531 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2532 fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2533 fHistograms->FillHistogram("ESD_Background_Pt", momentumVectorbackgroundCandidate.Pt());
2534 fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2535 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2536 fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2537 fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2538 fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2539 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2540 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2541 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2542 fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
111d75df 2543 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
6de3471d 2544
2545
2546 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(currentEventGoodV02.GetEta())<0.9 ){
2547 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
2548 fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
2549 }
2550
2551 fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
2552 fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
2553 fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin), momentumVectorbackgroundCandidate.Pt());
2554 fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
2555 fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
2556 fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
2557 fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
2558 fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2559 fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2560 fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
2561 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2562 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
2563
2564 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(currentEventGoodV02.GetEta())<0.9 ){
2565 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2566 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
2567 }
111d75df 2568 }
3c45d101 2569 if(alfa<0.1){
2570 fHistograms->FillHistogram("ESD_Background_InvMass_vs_E_alpha",massBG ,backgroundCandidate->GetE());
2571 }
2572
2e2da371 2573 }
111d75df 2574 //}
5ce758b0 2575 delete backgroundCandidate;
2576 }
2577 }
2578 }
111d75df 2579 delete random;
5ce758b0 2580 }
111d75df 2581 else{ // means no rotation
2582 AliGammaConversionBGHandler::GammaConversionVertex *bgEventVertex = NULL;
2583
2584 if(fUseTrackMultiplicityForBG){
2585 // cout<<"Using charged track multiplicity for background calculation"<<endl;
2586 for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
5ce758b0 2587
111d75df 2588 AliGammaConversionKFVector * previousEventV0s = bgHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);//fV0Reader->GetBGGoodV0s(nEventsInBG);
2589
5ce758b0 2590 if(fMoveParticleAccordingToVertex == kTRUE){
2591 bgEventVertex = bgHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
2592 }
2593
2594 for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2595 AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
2596 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2597 AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
111d75df 2598 AliKFParticle previousGoodV0test = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
5ce758b0 2599
111d75df 2600 //cout<<"Primary Vertex event: ["<<fESDEvent->GetPrimaryVertex()->GetX()<<","<<fESDEvent->GetPrimaryVertex()->GetY()<<","<<fESDEvent->GetPrimaryVertex()->GetZ()<<"]"<<endl;
2601 //cout<<"BG prim Vertex event: ["<<bgEventVertex->fX<<","<<bgEventVertex->fY<<","<<bgEventVertex->fZ<<"]"<<endl;
2602
2603 //cout<<"XYZ of particle before transport: ["<<previousGoodV0.X()<<","<<previousGoodV0.Y()<<","<<previousGoodV0.Z()<<"]"<<endl;
5ce758b0 2604 if(fMoveParticleAccordingToVertex == kTRUE){
2605 MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
2606 }
111d75df 2607 //cout<<"XYZ of particle after transport: ["<<previousGoodV0.X()<<","<<previousGoodV0.Y()<<","<<previousGoodV0.Z()<<"]"<<endl;
5ce758b0 2608
2609 AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
111d75df 2610
5ce758b0 2611 Double_t massBG =0.;
2612 Double_t widthBG = 0.;
2613 Double_t chi2BG =10000.;
2614 backgroundCandidate->GetMass(massBG,widthBG);
111d75df 2615 // if(backgroundCandidate->GetNDF()>0){
2616 // chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
2617 chi2BG = backgroundCandidate->GetChi2();
2618 if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
2619
2620 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2621 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2622
2623 Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2624
2625 Double_t rapidity;
2626
2627 if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() <= 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() <= 0){
2628 cout << "Error: |Pz| > E !!!! " << endl;
2629 rapidity=0;
2630 } else {
2631 rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2632 }
dc2883e4 2633 if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ){
2634 delete backgroundCandidate;
2635 continue; // rapidity cut
2636 }
2637
2638
111d75df 2639 Double_t alfa=0.0;
2640 if( (currentEventGoodV0.GetE()+previousGoodV0.GetE()) != 0){
2641 alfa=TMath::Abs((currentEventGoodV0.GetE()-previousGoodV0.GetE())
2642 /(currentEventGoodV0.GetE()+previousGoodV0.GetE()));
2643 }
2644
2645
2646 if(openingAngleBG < fMinOpeningAngleGhostCut ){
2647 delete backgroundCandidate;
2648 continue; // minimum opening angle to avoid using ghosttracks
2649 }
2650
2651 // original
111d75df 2652 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
6de3471d 2653 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2654 fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2655 fHistograms->FillHistogram("ESD_Background_Pt", momentumVectorbackgroundCandidate.Pt());
2656 fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2657 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2658 fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2659 fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2660 fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2661 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2662 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2663 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2664 fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
111d75df 2665 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
6de3471d 2666
2667
2668 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2669 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
2670 fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
2671 }
2672
2673 // test
2674 fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
2675 fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
2676 fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin), momentumVectorbackgroundCandidate.Pt());
2677 fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
2678 fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
2679 fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
2680 fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
2681 fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2682 fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2683 fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
2684 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2685 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
2686
2687 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2688 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2689 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
2690 }
2691 // }
111d75df 2692 }
3c45d101 2693 if(alfa<0.1){
2694 fHistograms->FillHistogram("ESD_Background_InvMass_vs_E_alpha",massBG ,backgroundCandidate->GetE());
2695 }
2696
111d75df 2697 }
2698 delete backgroundCandidate;
2699 }
2700 }
2701 }
2702 }
2703 else{ // means using #V0s for multiplicity
2704
2705 // cout<<"Using the v0 multiplicity to calculate background"<<endl;
2706
2707 fHistograms->FillHistogram("ESD_Background_z_m",zbin,mbin);
2708 fHistograms->FillHistogram("ESD_Mother_multpilicityVSv0s",fV0Reader->CountESDTracks(),fV0Reader->GetNumberOfV0s());
2709
2710 for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2711 AliGammaConversionKFVector * previousEventV0s = bgHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);// fV0Reader->GetBGGoodV0s(nEventsInBG);
2712 if(previousEventV0s){
2713
2714 if(fMoveParticleAccordingToVertex == kTRUE){
2715 bgEventVertex = bgHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
2716 }
2717
2718 for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2719 AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
2720 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2721 AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2722
2723 if(fMoveParticleAccordingToVertex == kTRUE){
2724 MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
5ce758b0 2725 }
111d75df 2726
2727 AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
2728 Double_t massBG =0.;
2729 Double_t widthBG = 0.;
2730 Double_t chi2BG =10000.;
2731 backgroundCandidate->GetMass(massBG,widthBG);
2732 /* if(backgroundCandidate->GetNDF()>0){
2733 chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
2734 {//remember to remove
2735 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2736 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2737
2738 Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2739 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle_nochi2", openingAngleBG);
2740 }
2741 */
2e2da371 2742 chi2BG = backgroundCandidate->GetChi2();
2743 if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
5ce758b0 2744 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2745 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2746
2747 Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2748
2749 Double_t rapidity;
2750 if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
2751 else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2752
dc2883e4 2753 if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ){
2754 delete backgroundCandidate;
2755 continue; // rapidity cut
2756 }
2757
2758
5ce758b0 2759 Double_t alfa=0.0;
2760 if( (currentEventGoodV0.GetE()+previousGoodV0.GetE()) != 0){
2761 alfa=TMath::Abs((currentEventGoodV0.GetE()-previousGoodV0.GetE())
2762 /(currentEventGoodV0.GetE()+previousGoodV0.GetE()));
2763 }
2764
2765
2766 if(openingAngleBG < fMinOpeningAngleGhostCut ){
2767 delete backgroundCandidate;
2768 continue; // minimum opening angle to avoid using ghosttracks
2769 }
2770
5ce758b0 2771 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
6de3471d 2772 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2773 fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2774 fHistograms->FillHistogram("ESD_Background_Pt", momentumVectorbackgroundCandidate.Pt());
2775 fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2776 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2777 fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2778 fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2779 fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2780 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2781 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2782 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2783 fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
2784
2785
5ce758b0 2786 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
3c45d101 2787
6de3471d 2788
2789 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2790 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
2791 fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
2792 }
2793
2794 if(massBG>0.5 && massBG<0.6){
2795 fHistograms->FillHistogram("ESD_Background_alfa_pt0506",momentumVectorbackgroundCandidate.Pt(),alfa);
2796 }
2797 if(massBG>0.3 && massBG<0.4){
2798 fHistograms->FillHistogram("ESD_Background_alfa_pt0304",momentumVectorbackgroundCandidate.Pt(),alfa);
2799 }
2800
2801 // test
2802 fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
2803 fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
2804 fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin), momentumVectorbackgroundCandidate.Pt());
2805 fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
2806 fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
2807 fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
2808 fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
2809 fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2810 fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2811 fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
2812 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2813 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
2814
2815 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2816 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2817 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
2818 }
5ce758b0 2819 }
5ce758b0 2820
6de3471d 2821 if(alfa<0.1){
2822 fHistograms->FillHistogram("ESD_Background_InvMass_vs_E_alpha",massBG ,backgroundCandidate->GetE());
5ce758b0 2823 }
2e2da371 2824 // }
5ce758b0 2825 }
111d75df 2826 delete backgroundCandidate;
2827 }
5ce758b0 2828 }
2829 }
d7d7e825 2830 }
111d75df 2831 } // end else (means use #v0s as multiplicity)
2832 } // end no rotation
d7d7e825 2833}
2834
2835
d7d7e825 2836void AliAnalysisTaskGammaConversion::ProcessGammasForGammaJetAnalysis(){
2837 //ProcessGammasForGammaJetAnalysis
a0b94e5c 2838
d7d7e825 2839 Double_t distIsoMin;
a0b94e5c 2840
d7d7e825 2841 CreateListOfChargedParticles();
a0b94e5c 2842
2843
6c84d371 2844 // for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
2845 for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
2846 AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
01b7fdcc 2847 TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
a0b94e5c 2848
01b7fdcc 2849 if( momentumVectorCurrentGamma.Pt()> fMinPtForGammaJet){
d7d7e825 2850 distIsoMin=GetMinimumDistanceToCharge(gammaIndex);
a0b94e5c 2851
d7d7e825 2852 if (distIsoMin > fMinIsoConeSize && fLeadingChargedIndex>=0){
01b7fdcc 2853 CalculateJetCone(gammaIndex);
d7d7e825 2854 }
2855 }
2856 }
2857}
2858
6272370b 2859//____________________________________________________________________
2860Bool_t AliAnalysisTaskGammaConversion::IsGoodImpPar(AliESDtrack *const track)
2861{
2862//
2863// check whether particle has good DCAr(Pt) impact
2864// parameter. Only for TPC+ITS tracks (7*sigma cut)
2865// Origin: Andrea Dainese
2866//
2867
2868Float_t d0z0[2],covd0z0[3];
2869track->GetImpactParameters(d0z0,covd0z0);
2870Float_t sigma= 0.0050+0.0060/TMath::Power(track->Pt(),0.9);
2871Float_t d0max = 7.*sigma;
2872if(TMath::Abs(d0z0[0]) < d0max) return kTRUE;
2873
2874return kFALSE;
2875}
2876
2877
d7d7e825 2878void AliAnalysisTaskGammaConversion::CreateListOfChargedParticles(){
2879 // CreateListOfChargedParticles
a0b94e5c 2880
d7d7e825 2881 fESDEvent = fV0Reader->GetESDEvent();
037dc2db 2882 Int_t numberOfESDTracks=0;
d7d7e825 2883 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
2884 AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
a0b94e5c 2885
d7d7e825 2886 if(!curTrack){
2887 continue;
2888 }
d707e3cf 2889 // Not needed if Standard function used.
2890// if(!IsGoodImpPar(curTrack)){
2891// continue;
2892// }
a0b94e5c 2893
d7d7e825 2894 if(fEsdTrackCuts->AcceptTrack(curTrack) ){
6c84d371 2895 new((*fChargedParticles)[fChargedParticles->GetEntriesFast()]) AliESDtrack(*curTrack);
2896 // fChargedParticles.push_back(curTrack);
d7d7e825 2897 fChargedParticlesId.push_back(iTracks);
037dc2db 2898 numberOfESDTracks++;
d7d7e825 2899 }
2900 }
037dc2db 2901 fHistograms->FillHistogram("ESD_NumberOfGoodESDTracks",numberOfESDTracks);
e40fd7e2 2902
2903 if (fV0Reader->GetNumberOfContributorsVtx()>=1){
2904 fHistograms->FillHistogram("ESD_NumberOfGoodESDTracksVtx",numberOfESDTracks);
2905 }