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