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