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