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