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