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