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