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