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