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