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