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