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