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