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