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