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