]> 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
Remove unnecessary print from QA class
[u/mrichter/AliRoot.git] / PWG4 / GammaConv / AliAnalysisTaskGammaConversion.cxx
CommitLineData
d7d7e825 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: Ana Marin, Kathrin Koch, Kenneth Aamodt *
5 * Version 1.1 *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16////////////////////////////////////////////////
17//---------------------------------------------
18// Class used to do analysis on conversion pairs
19//---------------------------------------------
20////////////////////////////////////////////////
21
22// root
23#include <TChain.h>
24
25// analysis
26#include "AliAnalysisTaskGammaConversion.h"
27#include "AliStack.h"
28#include "AliLog.h"
29#include "AliESDtrackCuts.h"
30#include "TNtuple.h"
4a6157dc 31//#include "AliCFManager.h" // for CF
32//#include "AliCFContainer.h" // for CF
9c1cb6f7 33#include "AliESDInputHandler.h"
34#include "AliAnalysisManager.h"
6441e967 35#include "AliAODPWG4Particle.h"
36#include "AliAODPWG4ParticleCorrelation.h"
4a6157dc 37#include "AliGammaConversionAODObject.h"
6441e967 38#include "AliAODConversionParticle.h"
5e55d806 39#include "AliGammaConversionBGHandler.h"
6272370b 40#include "AliESDCaloCluster.h" // for combining PHOS and GammaConv
9c1cb6f7 41#include "AliKFVertex.h"
3c45d101 42#include "AliGenPythiaEventHeader.h"
43#include "AliGenDPMjetEventHeader.h"
44#include "AliGenEventHeader.h"
e60f3265 45#include <AliMCEventHandler.h>
77ac6f3e 46#include "TRandom3.h"
18a1e325 47#include "AliTriggerAnalysis.h"
b6dd6ad2 48#include "AliCentrality.h"
77ac6f3e 49
d707e3cf 50class AliESDTrackCuts;
4a6157dc 51class AliCFContainer;
52class AliCFManager;
d7d7e825 53class AliKFVertex;
54class AliAODHandler;
55class AliAODEvent;
56class ALiESDEvent;
57class AliMCEvent;
58class AliMCEventHandler;
59class AliESDInputHandler;
60class AliAnalysisManager;
61class Riostream;
62class TFile;
63class TInterpreter;
64class TSystem;
65class TROOT;
66
67ClassImp(AliAnalysisTaskGammaConversion)
68
69
70AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion():
71AliAnalysisTaskSE(),
72 fV0Reader(NULL),
73 fStack(NULL),
a0b94e5c 74 fMCTruth(NULL), // for CF
4a6157dc 75 fGCMCEvent(NULL), // for CF
d7d7e825 76 fESDEvent(NULL),
77 fOutputContainer(NULL),
a0b94e5c 78 fCFManager(0x0), // for CF
d7d7e825 79 fHistograms(NULL),
b5832f95 80 fTriggerCINT1B(kFALSE),
d7d7e825 81 fDoMCTruth(kFALSE),
82 fDoNeutralMeson(kFALSE),
6272370b 83 fDoOmegaMeson(kFALSE),
d7d7e825 84 fDoJet(kFALSE),
85 fDoChic(kFALSE),
9c1cb6f7 86 fRecalculateV0ForGamma(kFALSE),
6c84d371 87 fKFReconstructedGammasTClone(NULL),
6272370b 88 fKFReconstructedPi0sTClone(NULL),
9c1cb6f7 89 fKFRecalculatedGammasTClone(NULL),
6c84d371 90 fCurrentEventPosElectronTClone(NULL),
91 fCurrentEventNegElectronTClone(NULL),
92 fKFReconstructedGammasCutTClone(NULL),
93 fPreviousEventTLVNegElectronTClone(NULL),
94 fPreviousEventTLVPosElectronTClone(NULL),
d7d7e825 95 fElectronv1(),
96 fElectronv2(),
6272370b 97 fGammav1(),
98 fGammav2(),
9c1cb6f7 99 fElectronRecalculatedv1(),
100 fElectronRecalculatedv2(),
d7d7e825 101 fElectronMass(-1),
102 fGammaMass(-1),
103 fPi0Mass(-1),
104 fEtaMass(-1),
105 fGammaWidth(-1),
106 fPi0Width(-1),
107 fEtaWidth(-1),
108 fMinOpeningAngleGhostCut(0.),
9640a3d1 109 fEsdTrackCuts(NULL),
d7d7e825 110 fCalculateBackground(kFALSE),
111 fWriteNtuple(kFALSE),
112 fGammaNtuple(NULL),
113 fNeutralMesonNtuple(NULL),
114 fTotalNumberOfAddedNtupleEntries(0),
6c84d371 115 fChargedParticles(NULL),
d7d7e825 116 fChargedParticlesId(),
117 fGammaPtHighest(0.),
118 fMinPtForGammaJet(1.),
119 fMinIsoConeSize(0.2),
120 fMinPtIsoCone(0.7),
121 fMinPtGamChargedCorr(0.5),
122 fMinPtJetCone(0.5),
123 fLeadingChargedIndex(-1),
9640a3d1 124 fLowPtMapping(1.),
125 fHighPtMapping(3.),
1e7846f4 126 fDoCF(kFALSE),
04bf4381 127 fAODGamma(NULL),
128 fAODPi0(NULL),
129 fAODOmega(NULL),
037dc2db 130 fAODBranchName("GammaConv"),
6441e967 131 fOutputAODClassName("AliAODConversionParticle"),
63477d82 132 fKFCreateAOD(kTRUE),
d765d400 133 fKFForceAOD(kFALSE),
332f1f44 134 fKFDeltaAODFileName(""),
037dc2db 135 fDoNeutralMesonV0MCCheck(kFALSE),
5ce758b0 136 fUseTrackMultiplicityForBG(kTRUE),
137 fMoveParticleAccordingToVertex(kFALSE),
2e2da371 138 fApplyChi2Cut(kFALSE),
77ac6f3e 139 fNRandomEventsForBG(15),
140 fNDegreesPMBackground(15),
111d75df 141 fDoRotation(kTRUE),
142 fCheckBGProbability(kTRUE),
2f426a53 143 fKFReconstructedGammasV0Index(),
18a1e325 144 fRemovePileUp(kFALSE),
145 fSelectV0AND(kFALSE),
6746e1e1 146 fTriggerAnalysis(NULL),
147 fMultiplicity(0),
148 fUseMultiplicity(0),
64f4118c 149 fUseMultiplicityBin(0),
150 fUseCentrality(0),
151 fUseCentralityBin(0)
d7d7e825 152{
153 // Default constructor
5a34881d 154
155 /* Kenneth: the default constructor should not have any define input/output or the call to SetESDtrackCuts
d7d7e825 156 // Common I/O in slot 0
157 DefineInput (0, TChain::Class());
158 DefineOutput(0, TTree::Class());
159
160 // Your private output
161 DefineOutput(1, TList::Class());
a0b94e5c 162
d7d7e825 163 // Define standard ESD track cuts for Gamma-hadron correlation
164 SetESDtrackCuts();
5a34881d 165 */
d7d7e825 166}
167
168AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):
169 AliAnalysisTaskSE(name),
170 fV0Reader(NULL),
171 fStack(NULL),
a0b94e5c 172 fMCTruth(NULL), // for CF
4a6157dc 173 fGCMCEvent(NULL), // for CF
d7d7e825 174 fESDEvent(NULL),
175 fOutputContainer(0x0),
a0b94e5c 176 fCFManager(0x0), // for CF
d7d7e825 177 fHistograms(NULL),
b5832f95 178 fTriggerCINT1B(kFALSE),
d7d7e825 179 fDoMCTruth(kFALSE),
180 fDoNeutralMeson(kFALSE),
6272370b 181 fDoOmegaMeson(kFALSE),
d7d7e825 182 fDoJet(kFALSE),
183 fDoChic(kFALSE),
9c1cb6f7 184 fRecalculateV0ForGamma(kFALSE),
6c84d371 185 fKFReconstructedGammasTClone(NULL),
6272370b 186 fKFReconstructedPi0sTClone(NULL),
9c1cb6f7 187 fKFRecalculatedGammasTClone(NULL),
6c84d371 188 fCurrentEventPosElectronTClone(NULL),
189 fCurrentEventNegElectronTClone(NULL),
190 fKFReconstructedGammasCutTClone(NULL),
191 fPreviousEventTLVNegElectronTClone(NULL),
192 fPreviousEventTLVPosElectronTClone(NULL),
d7d7e825 193 fElectronv1(),
194 fElectronv2(),
6272370b 195 fGammav1(),
196 fGammav2(),
9c1cb6f7 197 fElectronRecalculatedv1(),
198 fElectronRecalculatedv2(),
d7d7e825 199 fElectronMass(-1),
200 fGammaMass(-1),
201 fPi0Mass(-1),
202 fEtaMass(-1),
203 fGammaWidth(-1),
204 fPi0Width(-1),
205 fEtaWidth(-1),
206 fMinOpeningAngleGhostCut(0.),
9640a3d1 207 fEsdTrackCuts(NULL),
d7d7e825 208 fCalculateBackground(kFALSE),
209 fWriteNtuple(kFALSE),
210 fGammaNtuple(NULL),
211 fNeutralMesonNtuple(NULL),
212 fTotalNumberOfAddedNtupleEntries(0),
6c84d371 213 fChargedParticles(NULL),
d7d7e825 214 fChargedParticlesId(),
215 fGammaPtHighest(0.),
216 fMinPtForGammaJet(1.),
217 fMinIsoConeSize(0.2),
218 fMinPtIsoCone(0.7),
219 fMinPtGamChargedCorr(0.5),
220 fMinPtJetCone(0.5),
221 fLeadingChargedIndex(-1),
9640a3d1 222 fLowPtMapping(1.),
223 fHighPtMapping(3.),
c59360eb 224 fDoCF(kFALSE),
225 fAODGamma(NULL),
226 fAODPi0(NULL),
227 fAODOmega(NULL),
228 fAODBranchName("GammaConv"),
6441e967 229 fOutputAODClassName("AliAODConversionParticle"),
63477d82 230 fKFCreateAOD(kTRUE),
c59360eb 231 fKFForceAOD(kFALSE),
232 fKFDeltaAODFileName(""),
233 fDoNeutralMesonV0MCCheck(kFALSE),
5ce758b0 234 fUseTrackMultiplicityForBG(kTRUE),
235 fMoveParticleAccordingToVertex(kFALSE),
2e2da371 236 fApplyChi2Cut(kFALSE),
77ac6f3e 237 fNRandomEventsForBG(15),
238 fNDegreesPMBackground(15),
111d75df 239 fDoRotation(kTRUE),
240 fCheckBGProbability(kTRUE),
2f426a53 241 fKFReconstructedGammasV0Index(),
18a1e325 242 fRemovePileUp(kFALSE),
243 fSelectV0AND(kFALSE),
6746e1e1 244 fTriggerAnalysis(NULL),
245 fMultiplicity(0),
246 fUseMultiplicity(0),
64f4118c 247 fUseMultiplicityBin(0),
248 fUseCentrality(0),
249 fUseCentralityBin(0)
d7d7e825 250{
251 // Common I/O in slot 0
252 DefineInput (0, TChain::Class());
253 DefineOutput(0, TTree::Class());
254
255 // Your private output
256 DefineOutput(1, TList::Class());
a0b94e5c 257 DefineOutput(2, AliCFContainer::Class()); // for CF
258
259
d7d7e825 260 // Define standard ESD track cuts for Gamma-hadron correlation
261 SetESDtrackCuts();
9c1cb6f7 262
d7d7e825 263}
264
265AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion()
266{
267 // Remove all pointers
268
269 if(fOutputContainer){
270 fOutputContainer->Clear() ;
271 delete fOutputContainer ;
272 }
273 if(fHistograms){
274 delete fHistograms;
275 }
276 if(fV0Reader){
277 delete fV0Reader;
278 }
a0b94e5c 279
280 // for CF
281 if(fCFManager){
282 delete fCFManager;
283 }
9640a3d1 284
285 if(fEsdTrackCuts){
286 delete fEsdTrackCuts;
287 }
288
04bf4381 289 //Delete AODs
290 if (fAODGamma) {
291 fAODGamma->Clear();
292 delete fAODGamma;
d7d7e825 293 }
04bf4381 294 fAODGamma = NULL;
295
296 if (fAODPi0) {
297 fAODPi0->Clear();
298 delete fAODPi0;
299 }
300 fAODPi0 = NULL;
d7d7e825 301
04bf4381 302 if (fAODOmega) {
303 fAODOmega->Clear();
304 delete fAODOmega;
305 }
306 fAODOmega = NULL;
307
18a1e325 308 if(fTriggerAnalysis) {
309 delete fTriggerAnalysis;
310 }
311
312
04bf4381 313}
d7d7e825 314
12464034 315
d7d7e825 316void AliAnalysisTaskGammaConversion::Init()
317{
318 // Initialization
4a6157dc 319 // AliLog::SetGlobalLogLevel(AliLog::kError);
d7d7e825 320}
321void AliAnalysisTaskGammaConversion::SetESDtrackCuts()
322{
323 // SetESDtrackCuts
9640a3d1 324 if (fEsdTrackCuts!=NULL){
325 delete fEsdTrackCuts;
326 }
d7d7e825 327 fEsdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
a0b94e5c 328 //standard cuts from:
329 //http://aliceinfo.cern.ch/alicvs/viewvc/PWG0/dNdEta/CreateCuts.C?revision=1.4&view=markup
ebcfaa7e 330
331 // Cuts used up to 3rd of March
332
333 // fEsdTrackCuts->SetMinNClustersTPC(50);
334// fEsdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
335// fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
336// fEsdTrackCuts->SetRequireITSRefit(kTRUE);
337// fEsdTrackCuts->SetMaxNsigmaToVertex(3);
338// fEsdTrackCuts->SetRequireSigmaToVertex(kTRUE);
339
340 //------- To be tested-----------
d707e3cf 341 // Cuts used up to 26th of Agost
342// Int_t minNClustersTPC = 70;
343// Double_t maxChi2PerClusterTPC = 4.0;
344// Double_t maxDCAtoVertexXY = 2.4; // cm
345// Double_t maxDCAtoVertexZ = 3.2; // cm
346// fEsdTrackCuts->SetRequireSigmaToVertex(kFALSE);
347// fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
348// fEsdTrackCuts->SetRequireITSRefit(kTRUE);
349// // fEsdTrackCuts->SetRequireTPCStandAlone(kTRUE);
350// fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
351// fEsdTrackCuts->SetMinNClustersTPC(minNClustersTPC);
352// fEsdTrackCuts->SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC);
353// fEsdTrackCuts->SetMaxDCAToVertexXY(maxDCAtoVertexXY);
354// fEsdTrackCuts->SetMaxDCAToVertexZ(maxDCAtoVertexZ);
355// fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
356// fEsdTrackCuts->SetEtaRange(-0.8, 0.8);
357// fEsdTrackCuts->SetPtRange(0.15);
358
359// fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst);
360
361
362// Using standard function for setting Cuts
363 Bool_t selectPrimaries=kTRUE;
406a77e0 364 fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
64f4118c 365 fEsdTrackCuts->SetMaxDCAToVertexZ(2);
d707e3cf 366 fEsdTrackCuts->SetEtaRange(-0.8, 0.8);
367 fEsdTrackCuts->SetPtRange(0.15);
ebcfaa7e 368
369 //----- From Jacek 10.03.03 ------------------/
370// minNClustersTPC = 70;
371// maxChi2PerClusterTPC = 4.0;
372// maxDCAtoVertexXY = 2.4; // cm
373// maxDCAtoVertexZ = 3.2; // cm
374
375// esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
376// esdTrackCuts->SetRequireTPCRefit(kFALSE);
377// esdTrackCuts->SetRequireTPCStandAlone(kTRUE);
378// esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
379// esdTrackCuts->SetMinNClustersTPC(minNClustersTPC);
380// esdTrackCuts->SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC);
381// esdTrackCuts->SetMaxDCAToVertexXY(maxDCAtoVertexXY);
382// esdTrackCuts->SetMaxDCAToVertexZ(maxDCAtoVertexZ);
383// esdTrackCuts->SetDCAToVertex2D(kTRUE);
384
385
386
d7d7e825 387 // fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
037dc2db 388 // fV0Reader->SetESDtrackCuts(fEsdTrackCuts);
d7d7e825 389}
390
c00009fb 391void AliAnalysisTaskGammaConversion::UserExec(Option_t */*option*/)
d7d7e825 392{
393 // Execute analysis for current event
6272370b 394
9c1cb6f7 395 // Load the esdpid from the esdhandler if exists (tender was applied) otherwise set the Bethe Bloch parameters
3c45d101 396 Int_t eventQuality=-1;
9c1cb6f7 397
398 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
399 AliESDInputHandler *esdHandler=0x0;
400 if ( (esdHandler=dynamic_cast<AliESDInputHandler*>(man->GetInputEventHandler())) && esdHandler->GetESDpid() ){
401 AliV0Reader::SetESDpid(esdHandler->GetESDpid());
402 } else {
403 //load esd pid bethe bloch parameters depending on the existance of the MC handler
404 // yes: MC parameters
405 // no: data parameters
406 if (!AliV0Reader::GetESDpid()){
407 if (fMCEvent ) {
408 AliV0Reader::InitESDpid();
409 } else {
410 AliV0Reader::InitESDpid(1);
411 }
412 }
413 }
414
d765d400 415 //Must set fForceAOD to true for the AOD to get filled. Should only be done when running independent chain / train.
416 if(fKFForceAOD) {
417 if (!AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()) AliFatal("Cannot run ESD filter without an output event handler");
418 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
419 }
48682642 420
421 if(fV0Reader == NULL){
422 // Write warning here cuts and so on are default if this ever happens
6272370b 423 }
48682642 424
e60f3265 425 if (fMCEvent ) {
426 // To avoid crashes due to unzip errors. Sometimes the trees are not there.
427
428 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
3c45d101 429 if (!mcHandler){
430 AliError("Could not retrive MC event handler!");
431 return;
432
433 eventQuality=0;
434 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
435 }
436 if (!mcHandler->InitOk() ){
437 eventQuality=0;
438 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
439 return;
440 }
441 if (!mcHandler->TreeK() ){
442 eventQuality=0;
443 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
444 return;
445 }
446 if (!mcHandler->TreeTR() ) {
447 eventQuality=0;
448 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
449 return;
450 }
e60f3265 451 }
452
48682642 453 fV0Reader->SetInputAndMCEvent(InputEvent(), MCEvent());
454
455 fV0Reader->Initialize();
456 fDoMCTruth = fV0Reader->GetDoMCTruth();
457
12464034 458 if(fAODGamma) fAODGamma->Delete();
459 if(fAODPi0) fAODPi0->Delete();
460 if(fAODOmega) fAODOmega->Delete();
a0b94e5c 461
6c84d371 462 if(fKFReconstructedGammasTClone == NULL){
463 fKFReconstructedGammasTClone = new TClonesArray("AliKFParticle",0);
464 }
465 if(fCurrentEventPosElectronTClone == NULL){
466 fCurrentEventPosElectronTClone = new TClonesArray("AliESDtrack",0);
467 }
468 if(fCurrentEventNegElectronTClone == NULL){
469 fCurrentEventNegElectronTClone = new TClonesArray("AliESDtrack",0);
470 }
471 if(fKFReconstructedGammasCutTClone == NULL){
472 fKFReconstructedGammasCutTClone = new TClonesArray("AliKFParticle",0);
473 }
474 if(fPreviousEventTLVNegElectronTClone == NULL){
475 fPreviousEventTLVNegElectronTClone = new TClonesArray("TLorentzVector",0);
476 }
477 if(fPreviousEventTLVPosElectronTClone == NULL){
478 fPreviousEventTLVPosElectronTClone = new TClonesArray("TLorentzVector",0);
479 }
480 if(fChargedParticles == NULL){
481 fChargedParticles = new TClonesArray("AliESDtrack",0);
482 }
6272370b 483
484 if(fKFReconstructedPi0sTClone == NULL){
485 fKFReconstructedPi0sTClone = new TClonesArray("AliKFParticle",0);
486 }
9c1cb6f7 487
18a1e325 488 if(fKFRecalculatedGammasTClone == NULL){
9c1cb6f7 489 fKFRecalculatedGammasTClone = new TClonesArray("AliKFParticle",0);
490 }
491
18a1e325 492 if(fTriggerAnalysis== NULL){
493 fTriggerAnalysis = new AliTriggerAnalysis;
494 }
495
6c84d371 496 //clear TClones
1e7846f4 497 fKFReconstructedGammasTClone->Delete();
498 fCurrentEventPosElectronTClone->Delete();
499 fCurrentEventNegElectronTClone->Delete();
500 fKFReconstructedGammasCutTClone->Delete();
501 fPreviousEventTLVNegElectronTClone->Delete();
502 fPreviousEventTLVPosElectronTClone->Delete();
6272370b 503 fKFReconstructedPi0sTClone->Delete();
9c1cb6f7 504 fKFRecalculatedGammasTClone->Delete();
5e55d806 505
d7d7e825 506 //clear vectors
d7d7e825 507 fElectronv1.clear();
508 fElectronv2.clear();
6272370b 509
510 fGammav1.clear();
511 fGammav2.clear();
512
9c1cb6f7 513 fElectronRecalculatedv1.clear();
514 fElectronRecalculatedv2.clear();
515
d7d7e825 516
1e7846f4 517 fChargedParticles->Delete();
5e55d806 518
d7d7e825 519 fChargedParticlesId.clear();
037dc2db 520
521 fKFReconstructedGammasV0Index.clear();
a0b94e5c 522
d7d7e825 523 //Clear the data in the v0Reader
5e55d806 524 // fV0Reader->UpdateEventByEventData();
b5832f95 525
526 //Take Only events with proper trigger
5e55d806 527 /*
b5832f95 528 if(fTriggerCINT1B){
529 if(!fV0Reader->GetESDEvent()->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) return;
530 }
5e55d806 531 */
10e3319b 532
c8206114 533 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
534 // cout<< "Event not taken"<< endl;
3c45d101 535 eventQuality=1;
536 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
537 if(fDoMCTruth){
538 CheckMesonProcessTypeEventQuality(eventQuality);
539 }
c8206114 540 return; // aborts if the primary vertex does not have contributors.
541 }
542
543
885114d1 544 if(!fV0Reader->CheckForPrimaryVertexZ() ){
3c45d101 545 eventQuality=2;
546 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
547 if(fDoMCTruth){
548 CheckMesonProcessTypeEventQuality(eventQuality);
549 }
885114d1 550 return;
551 }
2f426a53 552
553 if(fRemovePileUp && fV0Reader->GetESDEvent()->IsPileupFromSPD()) {
0a3d9f49 554 eventQuality=4;
45441b78 555 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
2f426a53 556 return;
557 }
18a1e325 558
18a1e325 559 Bool_t v0A = fTriggerAnalysis->IsOfflineTriggerFired(fV0Reader->GetESDEvent(), AliTriggerAnalysis::kV0A);
560 Bool_t v0C = fTriggerAnalysis->IsOfflineTriggerFired(fV0Reader->GetESDEvent(), AliTriggerAnalysis::kV0C);
561 Bool_t v0AND = v0A && v0C;
562
563 if(fSelectV0AND && !v0AND){
564 eventQuality=5;
45441b78 565 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
18a1e325 566 return;
567 }
6746e1e1 568 fMultiplicity = fEsdTrackCuts->CountAcceptedTracks(fV0Reader->GetESDEvent());
569
570 if( CalculateMultiplicityBin() != fUseMultiplicityBin){
571 eventQuality=6;
572 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
573 return;
574 }
18a1e325 575
64f4118c 576 if(fV0Reader->GetIsHeavyIon()){
577 if(fUseCentrality>0){
b6dd6ad2 578 AliCentrality *esdCentrality = fV0Reader->GetESDEvent()->GetCentrality();
64f4118c 579 Int_t centralityC = -1;
580
581 if(fUseCentrality==1){
582 centralityC = esdCentrality->GetCentralityClass10("V0M");
583 if( centralityC != fUseCentralityBin ){
584 eventQuality=7;
585 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
586 return;
587 }
588 }
589
590 if(fUseCentrality==2){
591 centralityC = esdCentrality->GetCentralityClass10("CL1");
592 if( centralityC != fUseCentralityBin ){
593 eventQuality=7;
594 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
595 return;
596 }
597 }
598 }
599 }
3c45d101 600 eventQuality=3;
601 fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
885114d1 602
6746e1e1 603
604
605 fHistograms->FillHistogram("ESD_NumberOfGoodESDTracks",fMultiplicity);
606 if (fV0Reader->GetNumberOfContributorsVtx()>=1){
607 fHistograms->FillHistogram("ESD_NumberOfGoodESDTracksVtx",fMultiplicity);
608 }
609
610
611
d7d7e825 612 // Process the MC information
613 if(fDoMCTruth){
614 ProcessMCData();
615 }
a0b94e5c 616
d7d7e825 617 //Process the v0 information with no cuts
618 ProcessV0sNoCut();
cb90a330 619
d7d7e825 620 // Process the v0 information
621 ProcessV0s();
cb90a330 622
5e55d806 623
d7d7e825 624 //Fill Gamma AOD
63477d82 625 if(fKFCreateAOD) {
626 FillAODWithConversionGammas() ;
627 }
628
629
d7d7e825 630 // Process reconstructed gammas
631 if(fDoNeutralMeson == kTRUE){
a0b94e5c 632 ProcessGammasForNeutralMesonAnalysis();
77880bd8 633
48682642 634 }
635
77880bd8 636 if(fDoMCTruth == kTRUE){
637 CheckV0Efficiency();
638 }
d7d7e825 639 //Process reconstructed gammas electrons for Chi_c Analysis
6c84d371 640 if(fDoChic == kTRUE){
d7d7e825 641 ProcessGammaElectronsForChicAnalysis();
642 }
643 // Process reconstructed gammas for gamma Jet/hadron correlations
644 if(fDoJet == kTRUE){
645 ProcessGammasForGammaJetAnalysis();
646 }
a0b94e5c 647
5e55d806 648 //calculate background if flag is set
649 if(fCalculateBackground){
650 CalculateBackground();
651 }
48682642 652
653 if(fDoNeutralMeson == kTRUE){
654// ProcessConvPHOSGammasForNeutralMesonAnalysis();
655 if(fDoOmegaMeson == kTRUE){
656 ProcessGammasForOmegaMesonAnalysis();
657 }
658 }
659
5e55d806 660 //Clear the data in the v0Reader
661 fV0Reader->UpdateEventByEventData();
9c1cb6f7 662 if(fRecalculateV0ForGamma==kTRUE){
663 RecalculateV0ForGamma();
664 }
d7d7e825 665 PostData(1, fOutputContainer);
a0b94e5c 666 PostData(2, fCFManager->GetParticleContainer()); // for CF
d7d7e825 667
668}
669
48682642 670// void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t *option){
671// // see header file for documentation
672// // printf(" ConnectInputData %s\n", GetName());
8a685cf3 673
48682642 674// AliAnalysisTaskSE::ConnectInputData(option);
8a685cf3 675
48682642 676// if(fV0Reader == NULL){
677// // Write warning here cuts and so on are default if this ever happens
678// }
679// fV0Reader->Initialize();
680// fDoMCTruth = fV0Reader->GetDoMCTruth();
681// }
d7d7e825 682
3c45d101 683void AliAnalysisTaskGammaConversion::CheckMesonProcessTypeEventQuality(Int_t evtQ){
77ac6f3e 684 // Check meson process type event quality
3c45d101 685 fStack= MCEvent()->Stack();
686 fGCMCEvent=MCEvent();
687
688 for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {
689 TParticle* particle = (TParticle *)fStack->Particle(iTracks);
690 if (!particle) {
691 //print warning here
692 continue;
693 }
694 if(particle->GetPdgCode()!=111){ //Pi0
695 continue;
696 }
697 if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() ) continue;
698 if(evtQ==1){
699 switch(GetProcessType(fGCMCEvent)){
700 case kProcSD:
701 fHistograms->FillHistogram("MC_SD_EvtQ1_Pi0_Pt", particle->Pt());
702 break;
703 case kProcDD:
704 fHistograms->FillHistogram("MC_DD_EvtQ1_Pi0_Pt", particle->Pt());
705 break;
706 case kProcND:
707 fHistograms->FillHistogram("MC_ND_EvtQ1_Pi0_Pt", particle->Pt());
708 break;
709 default:
710 AliError("Unknown Process");
711 }
712 }
713 if(evtQ==2){
714 switch(GetProcessType(fGCMCEvent)){
715 case kProcSD:
716 fHistograms->FillHistogram("MC_SD_EvtQ2_Pi0_Pt", particle->Pt());
717 break;
718 case kProcDD:
719 fHistograms->FillHistogram("MC_DD_EvtQ2_Pi0_Pt", particle->Pt());
720 break;
721 case kProcND:
722 fHistograms->FillHistogram("MC_ND_EvtQ2_Pi0_Pt", particle->Pt());
723 break;
724 default:
725 AliError("Unknown Process");
726 }
727 }
d7d7e825 728
729
3c45d101 730 }
731
732}
733
d7d7e825 734void AliAnalysisTaskGammaConversion::ProcessMCData(){
735 // see header file for documentation
48682642 736 //InputEvent(), MCEvent());
737 /* TestAnaMarin
d7d7e825 738 fStack = fV0Reader->GetMCStack();
a0b94e5c 739 fMCTruth = fV0Reader->GetMCTruth(); // for CF
4a6157dc 740 fGCMCEvent = fV0Reader->GetMCEvent(); // for CF
48682642 741 */
742 fStack= MCEvent()->Stack();
743 fGCMCEvent=MCEvent();
a0b94e5c 744
745 // for CF
1e7846f4 746 Double_t containerInput[3];
747 if(fDoCF){
748 if(!fGCMCEvent) cout << "NO MC INFO FOUND" << endl;
749 fCFManager->SetEventInfo(fGCMCEvent);
750 }
a0b94e5c 751 // end for CF
752
753
d7d7e825 754 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
755 return; // aborts if the primary vertex does not have contributors.
756 }
2f426a53 757 for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {
758 // for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {
d7d7e825 759 TParticle* particle = (TParticle *)fStack->Particle(iTracks);
bd6d9fa3 760
111d75df 761
bd6d9fa3 762
d7d7e825 763 if (!particle) {
764 //print warning here
765 continue;
766 }
a0b94e5c 767
d7d7e825 768 ///////////////////////Begin Chic Analysis/////////////////////////////
91b98865 769 if(fDoChic) {
770 if(particle->GetPdgCode() == 443){//Is JPsi
771 if(particle->GetNDaughters()==2){
772 if(TMath::Abs(fStack->Particle(particle->GetFirstDaughter())->GetPdgCode()) == 11 &&
773 TMath::Abs(fStack->Particle(particle->GetLastDaughter())->GetPdgCode()) == 11){
774
775 TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
776 TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
777 if(TMath::Abs(daug0->Eta()) < 0.9 && TMath::Abs(daug1->Eta()) < 0.9)
778 fHistograms->FillTable("Table_Electrons",3);//e+ e- from J/Psi inside acceptance
a0b94e5c 779
91b98865 780 if( TMath::Abs(daug0->Eta()) < 0.9){
781 if(daug0->GetPdgCode() == -11)
782 fHistograms->FillTable("Table_Electrons",1);//e+ from J/Psi inside acceptance
783 else
784 fHistograms->FillTable("Table_Electrons",2);//e- from J/Psi inside acceptance
a0b94e5c 785
91b98865 786 }
787 if(TMath::Abs(daug1->Eta()) < 0.9){
788 if(daug1->GetPdgCode() == -11)
789 fHistograms->FillTable("Table_Electrons",1);//e+ from J/Psi inside acceptance
790 else
791 fHistograms->FillTable("Table_Electrons",2);//e- from J/Psi inside acceptance
792 }
d7d7e825 793 }
794 }
795 }
91b98865 796 // const int CHI_C0 = 10441;
797 // const int CHI_C1 = 20443;
798 // const int CHI_C2 = 445
799 if(particle->GetPdgCode() == 22){//gamma from JPsi
800 if(particle->GetMother(0) > -1){
801 if(fStack->Particle(particle->GetMother(0))->GetPdgCode() == 10441 ||
802 fStack->Particle(particle->GetMother(0))->GetPdgCode() == 20443 ||
803 fStack->Particle(particle->GetMother(0))->GetPdgCode() == 445){
804 if(TMath::Abs(particle->Eta()) < 1.2)
805 fHistograms->FillTable("Table_Electrons",17);// gamma from chic inside accptance
806 }
d7d7e825 807 }
808 }
91b98865 809 if(particle->GetPdgCode() == 10441 || particle->GetPdgCode() == 20443 || particle->GetPdgCode() == 445){
810 if( particle->GetNDaughters() == 2){
811 TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
812 TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
a0b94e5c 813
91b98865 814 if( (daug0->GetPdgCode() == 443 || daug0->GetPdgCode() == 22) && (daug1->GetPdgCode() == 443 || daug1->GetPdgCode() == 22) ){
815 if( daug0->GetPdgCode() == 443){
816 TParticle* daugE0 = fStack->Particle(daug0->GetFirstDaughter());
817 TParticle* daugE1 = fStack->Particle(daug0->GetLastDaughter());
818 if( TMath::Abs(daug1->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
819 fHistograms->FillTable("Table_Electrons",18);
a0b94e5c 820
91b98865 821 }//if
822 else if (daug1->GetPdgCode() == 443){
823 TParticle* daugE0 = fStack->Particle(daug1->GetFirstDaughter());
824 TParticle* daugE1 = fStack->Particle(daug1->GetLastDaughter());
825 if( TMath::Abs(daug0->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
826 fHistograms->FillTable("Table_Electrons",18);
827 }//else if
828 }//gamma o Jpsi
829 }//GetNDaughters
830 }
d7d7e825 831 }
a0b94e5c 832
d7d7e825 833 /////////////////////End Chic Analysis////////////////////////////
a0b94e5c 834
835
dc2883e4 836 // if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() ) continue;
a0b94e5c 837
d7d7e825 838 if(particle->R()>fV0Reader->GetMaxRCut()) continue; // cuts on distance from collision point
839
840 Double_t tmpPhi=particle->Phi();
841
842 if(particle->Phi()> TMath::Pi()){
843 tmpPhi = particle->Phi()-(2*TMath::Pi());
844 }
845
846 Double_t rapidity;
847 if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){
848 rapidity=0;
849 }
850 else{
851 rapidity = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())));
852 }
853
111d75df 854
855
856 if(iTracks<=fStack->GetNprimary() ){
857 if ( particle->GetPdgCode()== -211 || particle->GetPdgCode()== 211 ||
858 particle->GetPdgCode()== 2212 || particle->GetPdgCode()==-2212 ||
859 particle->GetPdgCode()== 321 || particle->GetPdgCode()==-321 ){
dc2883e4 860 if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() ) continue;
3c45d101 861 fHistograms->FillHistogram("MC_PhysicalPrimaryCharged_Pt", particle->Pt());
111d75df 862 }
bd6d9fa3 863 }
111d75df 864
bd6d9fa3 865
866
d7d7e825 867 //process the gammas
868 if (particle->GetPdgCode() == 22){
dc2883e4 869 if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() ) continue;
a607c218 870
d7d7e825 871 if(particle->GetMother(0) >-1 && fStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){
872 continue; // no photon as mothers!
873 }
a607c218 874
d7d7e825 875 if(particle->GetMother(0) >= fStack->GetNprimary()){
876 continue; // the gamma has a mother, and it is not a primary particle
877 }
f8017b04 878
879 if(particle->GetMother(0) >-1){
880 fHistograms->FillHistogram("MC_DecayAllGamma_Pt", particle->Pt()); // All
881 switch(fStack->Particle(particle->GetMother(0))->GetPdgCode()){
882 case 111: // Pi0
883 fHistograms->FillHistogram("MC_DecayPi0Gamma_Pt", particle->Pt());
884 break;
885 case 113: // Rho0
886 fHistograms->FillHistogram("MC_DecayRho0Gamma_Pt", particle->Pt());
887 break;
888 case 221: // Eta
889 fHistograms->FillHistogram("MC_DecayEtaGamma_Pt", particle->Pt());
890 break;
891 case 223: // Omega
892 fHistograms->FillHistogram("MC_DecayOmegaGamma_Pt", particle->Pt());
893 break;
894 case 310: // K_s0
895 fHistograms->FillHistogram("MC_DecayK0sGamma_Pt", particle->Pt());
896 break;
897 case 331: // Eta'
898 fHistograms->FillHistogram("MC_DecayEtapGamma_Pt", particle->Pt());
899 break;
900 }
901 }
902
d7d7e825 903 fHistograms->FillHistogram("MC_allGamma_Energy", particle->Energy());
904 fHistograms->FillHistogram("MC_allGamma_Pt", particle->Pt());
905 fHistograms->FillHistogram("MC_allGamma_Eta", particle->Eta());
906 fHistograms->FillHistogram("MC_allGamma_Phi", tmpPhi);
907 fHistograms->FillHistogram("MC_allGamma_Rapid", rapidity);
908
a0b94e5c 909 // for CF
1e7846f4 910 if(fDoCF){
911 containerInput[0] = particle->Pt();
912 containerInput[1] = particle->Eta();
913 if(particle->GetMother(0) >=0){
914 containerInput[2] = fStack->Particle(particle->GetMother(0))->GetMass();
915 }
916 else{
917 containerInput[2]=-1;
918 }
a607c218 919
1e7846f4 920 fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated); // generated gamma
921 }
d7d7e825 922 if(particle->GetMother(0) < 0){ // direct gamma
923 fHistograms->FillHistogram("MC_allDirectGamma_Energy",particle->Energy());
924 fHistograms->FillHistogram("MC_allDirectGamma_Pt", particle->Pt());
925 fHistograms->FillHistogram("MC_allDirectGamma_Eta", particle->Eta());
926 fHistograms->FillHistogram("MC_allDirectGamma_Phi", tmpPhi);
927 fHistograms->FillHistogram("MC_allDirectGamma_Rapid", rapidity);
928 }
929
d7d7e825 930 // looking for conversion (electron + positron from pairbuilding (= 5) )
931 TParticle* ePos = NULL;
932 TParticle* eNeg = NULL;
933
934 if(particle->GetNDaughters() >= 2){
935 for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
936 TParticle *tmpDaughter = fStack->Particle(daughterIndex);
937 if(tmpDaughter->GetUniqueID() == 5){
938 if(tmpDaughter->GetPdgCode() == 11){
939 eNeg = tmpDaughter;
940 }
941 else if(tmpDaughter->GetPdgCode() == -11){
942 ePos = tmpDaughter;
943 }
944 }
945 }
946 }
947
948
949 if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
950 continue;
951 }
952
953
954 Double_t ePosPhi = ePos->Phi();
955 if(ePos->Phi()> TMath::Pi()) ePosPhi = ePos->Phi()-(2*TMath::Pi());
956
957 Double_t eNegPhi = eNeg->Phi();
958 if(eNeg->Phi()> TMath::Pi()) eNegPhi = eNeg->Phi()-(2*TMath::Pi());
959
960
961 if(ePos->Pt()<fV0Reader->GetPtCut() || eNeg->Pt()<fV0Reader->GetPtCut()){
962 continue; // no reconstruction below the Pt cut
963 }
a0b94e5c 964
d7d7e825 965 if(TMath::Abs(ePos->Eta())> fV0Reader->GetEtaCut() || TMath::Abs(eNeg->Eta())> fV0Reader->GetEtaCut()){
966 continue;
967 }
a0b94e5c 968
d7d7e825 969 if(ePos->R()>fV0Reader->GetMaxRCut()){
970 continue; // cuts on distance from collision point
971 }
a0b94e5c 972
973 if(TMath::Abs(ePos->Vz()) > fV0Reader->GetMaxZCut()){
d7d7e825 974 continue; // outside material
975 }
976
a0b94e5c 977
d7d7e825 978 if((TMath::Abs(ePos->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() > ePos->R()){
979 continue; // line cut to exclude regions where we do not reconstruct
980 }
981
a0b94e5c 982
983 // for CF
1e7846f4 984 if(fDoCF){
985 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructable); // reconstructable gamma
986 }
d7d7e825 987 fHistograms->FillHistogram("MC_ConvGamma_Energy", particle->Energy());
988 fHistograms->FillHistogram("MC_ConvGamma_Pt", particle->Pt());
989 fHistograms->FillHistogram("MC_ConvGamma_Eta", particle->Eta());
990 fHistograms->FillHistogram("MC_ConvGamma_Phi", tmpPhi);
991 fHistograms->FillHistogram("MC_ConvGamma_Rapid", rapidity);
992 fHistograms->FillHistogram("MC_ConvGamma_Pt_Eta", particle->Pt(),particle->Eta());
993
994 fHistograms->FillHistogram("MC_E_Energy", eNeg->Energy());
995 fHistograms->FillHistogram("MC_E_Pt", eNeg->Pt());
996 fHistograms->FillHistogram("MC_E_Eta", eNeg->Eta());
997 fHistograms->FillHistogram("MC_E_Phi", eNegPhi);
998
999 fHistograms->FillHistogram("MC_P_Energy", ePos->Energy());
1000 fHistograms->FillHistogram("MC_P_Pt", ePos->Pt());
1001 fHistograms->FillHistogram("MC_P_Eta", ePos->Eta());
1002 fHistograms->FillHistogram("MC_P_Phi", ePosPhi);
1003
1004
d7d7e825 1005 // begin Mapping
1006 Int_t rBin = fHistograms->GetRBin(ePos->R());
9640a3d1 1007 Int_t zBin = fHistograms->GetZBin(ePos->Vz());
d7d7e825 1008 Int_t phiBin = fHistograms->GetPhiBin(particle->Phi());
9c1cb6f7 1009 Double_t rFMD=30;
64f4118c 1010 Double_t rITSTPCMin=50;
1011 Double_t rITSTPCMax=80;
9640a3d1 1012
1013 TVector3 vtxPos(ePos->Vx(),ePos->Vy(),ePos->Vz());
1014
d7d7e825 1015 TString nameMCMappingPhiR="";
e158cbc3 1016 nameMCMappingPhiR.Form("MC_Conversion_Mapping_Phi%02d_R%02d",phiBin,rBin);
ebcfaa7e 1017 // fHistograms->FillHistogram(nameMCMappingPhiR, ePos->Vz(), particle->Eta());
d7d7e825 1018
1019 TString nameMCMappingPhi="";
e158cbc3 1020 nameMCMappingPhi.Form("MC_Conversion_Mapping_Phi%02d",phiBin);
9640a3d1 1021 // fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());
ebcfaa7e 1022 //fHistograms->FillHistogram(nameMCMappingPhi, ePos->Vz(), particle->Eta());
d7d7e825 1023
1024 TString nameMCMappingR="";
e158cbc3 1025 nameMCMappingR.Form("MC_Conversion_Mapping_R%02d",rBin);
9640a3d1 1026 // fHistograms->FillHistogram(nameMCMappingR, particle->Eta());
ebcfaa7e 1027 //fHistograms->FillHistogram(nameMCMappingR,ePos->Vz(), particle->Eta());
d7d7e825 1028
1029 TString nameMCMappingPhiInR="";
e158cbc3 1030 nameMCMappingPhiInR.Form("MC_Conversion_Mapping_Phi_in_R_%02d",rBin);
9640a3d1 1031 // fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
1032 fHistograms->FillHistogram(nameMCMappingPhiInR, vtxPos.Phi());
1033
1034 TString nameMCMappingZInR="";
1035 nameMCMappingZInR.Form("MC_Conversion_Mapping_Z_in_R_%02d",rBin);
1036 fHistograms->FillHistogram(nameMCMappingZInR,ePos->Vz() );
1037
1038
1039 TString nameMCMappingPhiInZ="";
1040 nameMCMappingPhiInZ.Form("MC_Conversion_Mapping_Phi_in_Z_%02d",zBin);
1041 // fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
1042 fHistograms->FillHistogram(nameMCMappingPhiInZ, vtxPos.Phi());
1043
9c1cb6f7 1044
1045 if(ePos->R()<rFMD){
1046 TString nameMCMappingFMDPhiInZ="";
1047 nameMCMappingFMDPhiInZ.Form("MC_Conversion_Mapping_FMD_Phi_in_Z_%02d",zBin);
1048 fHistograms->FillHistogram(nameMCMappingFMDPhiInZ, vtxPos.Phi());
1049 }
1050
64f4118c 1051 if(ePos->R()>rITSTPCMin && ePos->R()<rITSTPCMax){
1052 TString nameMCMappingITSTPCPhiInZ="";
1053 nameMCMappingITSTPCPhiInZ.Form("MC_Conversion_Mapping_ITSTPC_Phi_in_Z_%02d",zBin);
1054 fHistograms->FillHistogram(nameMCMappingITSTPCPhiInZ, vtxPos.Phi());
1055 }
1056
9640a3d1 1057 TString nameMCMappingRInZ="";
1058 nameMCMappingRInZ.Form("MC_Conversion_Mapping_R_in_Z_%02d",zBin);
1059 fHistograms->FillHistogram(nameMCMappingRInZ,ePos->R() );
1060
1061 if(particle->Pt() > fLowPtMapping && particle->Pt()< fHighPtMapping){
1062 TString nameMCMappingMidPtPhiInR="";
1063 nameMCMappingMidPtPhiInR.Form("MC_Conversion_Mapping_MidPt_Phi_in_R_%02d",rBin);
1064 fHistograms->FillHistogram(nameMCMappingMidPtPhiInR, vtxPos.Phi());
1065
1066 TString nameMCMappingMidPtZInR="";
1067 nameMCMappingMidPtZInR.Form("MC_Conversion_Mapping_MidPt_Z_in_R_%02d",rBin);
1068 fHistograms->FillHistogram(nameMCMappingMidPtZInR,ePos->Vz() );
1069
1070
1071 TString nameMCMappingMidPtPhiInZ="";
1072 nameMCMappingMidPtPhiInZ.Form("MC_Conversion_Mapping_MidPt_Phi_in_Z_%02d",zBin);
1073 fHistograms->FillHistogram(nameMCMappingMidPtPhiInZ, vtxPos.Phi());
9c1cb6f7 1074
1075
1076 if(ePos->R()<rFMD){
1077 TString nameMCMappingMidPtFMDPhiInZ="";
1078 nameMCMappingMidPtFMDPhiInZ.Form("MC_Conversion_Mapping_MidPt_FMD_Phi_in_Z_%02d",zBin);
1079 fHistograms->FillHistogram(nameMCMappingMidPtFMDPhiInZ, vtxPos.Phi());
1080 }
1081
9640a3d1 1082
1083 TString nameMCMappingMidPtRInZ="";
1084 nameMCMappingMidPtRInZ.Form("MC_Conversion_Mapping_MidPt_R_in_Z_%02d",zBin);
1085 fHistograms->FillHistogram(nameMCMappingMidPtRInZ,ePos->R() );
1086
1087 }
1088
d7d7e825 1089 //end mapping
1090
1091 fHistograms->FillHistogram("MC_Conversion_R",ePos->R());
1092 fHistograms->FillHistogram("MC_Conversion_ZR",ePos->Vz(),ePos->R());
1093 fHistograms->FillHistogram("MC_Conversion_XY",ePos->Vx(),ePos->Vy());
1094 fHistograms->FillHistogram("MC_Conversion_OpeningAngle",GetMCOpeningAngle(ePos, eNeg));
9640a3d1 1095 fHistograms->FillHistogram("MC_ConvGamma_E_AsymmetryP",particle->P(),eNeg->P()/particle->P());
1096 fHistograms->FillHistogram("MC_ConvGamma_P_AsymmetryP",particle->P(),ePos->P()/particle->P());
1097
1098
d7d7e825 1099 if(particle->GetMother(0) < 0){ // no mother = direct gamma, still inside converted
1100 fHistograms->FillHistogram("MC_ConvDirectGamma_Energy",particle->Energy());
1101 fHistograms->FillHistogram("MC_ConvDirectGamma_Pt", particle->Pt());
1102 fHistograms->FillHistogram("MC_ConvDirectGamma_Eta", particle->Eta());
1103 fHistograms->FillHistogram("MC_ConvDirectGamma_Phi", tmpPhi);
1104 fHistograms->FillHistogram("MC_ConvDirectGamma_Rapid", rapidity);
1105
1106 } // end direct gamma
1107 else{ // mother exits
6c84d371 1108 /* if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0
a0b94e5c 1109 fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S
1110 fStack->Particle(particle->GetMother(0))->GetPdgCode()==445 //chic2
1111 ){
1112 fMCGammaChic.push_back(particle);
1113 }
6c84d371 1114 */
d7d7e825 1115 } // end if mother exits
1116 } // end if particle is a photon
1117
1118
1119
1120 // process motherparticles (2 gammas as daughters)
1121 // the motherparticle had already to pass the R and the eta cut, but no line cut.
1122 // the line cut is just valid for the conversions!
1123
1124 if(particle->GetNDaughters() == 2){
1125
1126 TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());
1127 TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());
1128
1129 if(daughter0->GetPdgCode() != 22 || daughter1->GetPdgCode() != 22) continue; //check for gamma gamma daughters
dc2883e4 1130
1131 if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ) continue;
1132
d7d7e825 1133 // Check the acceptance for both gammas
01b7fdcc 1134 Bool_t gammaEtaCut = kTRUE;
1135 if(TMath::Abs(daughter0->Eta()) > fV0Reader->GetEtaCut() || TMath::Abs(daughter1->Eta()) > fV0Reader->GetEtaCut() ) gammaEtaCut = kFALSE;
1136 Bool_t gammaRCut = kTRUE;
1137 if(daughter0->R() > fV0Reader->GetMaxRCut() || daughter1->R() > fV0Reader->GetMaxRCut() ) gammaRCut = kFALSE;
d7d7e825 1138
1139
1140
1141 // check for conversions now -> have to pass eta, R and line cut!
1142 Bool_t daughter0Electron = kFALSE;
1143 Bool_t daughter0Positron = kFALSE;
1144 Bool_t daughter1Electron = kFALSE;
1145 Bool_t daughter1Positron = kFALSE;
1146
1147 if(daughter0->GetNDaughters() >= 2){ // first gamma
1148 for(Int_t TrackIndex=daughter0->GetFirstDaughter();TrackIndex<=daughter0->GetLastDaughter();TrackIndex++){
1149 TParticle *tmpDaughter = fStack->Particle(TrackIndex);
1150 if(tmpDaughter->GetUniqueID() == 5){
1151 if(tmpDaughter->GetPdgCode() == 11){
1152 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
1153 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
1154 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
1155 daughter0Electron = kTRUE;
1156 }
1157 }
1158
1159 }
1160 }
1161 else if(tmpDaughter->GetPdgCode() == -11){
1162 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
1163 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
1164 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
1165 daughter0Positron = kTRUE;
1166 }
a68437fb 1167 }
1168 }
d7d7e825 1169 }
1170 }
1171 }
1172 }
1173
1174
1175 if(daughter1->GetNDaughters() >= 2){ // second gamma
1176 for(Int_t TrackIndex=daughter1->GetFirstDaughter();TrackIndex<=daughter1->GetLastDaughter();TrackIndex++){
1177 TParticle *tmpDaughter = fStack->Particle(TrackIndex);
1178 if(tmpDaughter->GetUniqueID() == 5){
1179 if(tmpDaughter->GetPdgCode() == 11){
1180 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
1181 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
1182 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
1183 daughter1Electron = kTRUE;
1184 }
1185 }
1186
1187 }
1188 }
1189 else if(tmpDaughter->GetPdgCode() == -11){
1190 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
1191 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
1192 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
1193 daughter1Positron = kTRUE;
1194 }
1195 }
1196
1197 }
1198
1199 }
1200 }
1201 }
1202 }
1203
a0b94e5c 1204
d7d7e825 1205 if(particle->GetPdgCode()==111){ //Pi0
1206 if( iTracks >= fStack->GetNprimary()){
1207 fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());
1208 fHistograms->FillHistogram("MC_Pi0_Secondaries_Rapid", rapidity);
1209 fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);
1210 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());
1211 fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());
1212 fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());
1213 fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());
1214 fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1215 fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling
1216
01b7fdcc 1217 if(gammaEtaCut && gammaRCut){
d7d7e825 1218 //if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
1219 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1220 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
1221 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1222 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1223 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1224 }
1225 }
1226 }
1227 else{
1228 fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());
1229 fHistograms->FillHistogram("MC_Pi0_Rapid", rapidity);
1230 fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);
1231 fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());
dc2883e4 1232 fHistograms->FillHistogram("MC_Pi0_Pt_vs_Rapid", particle->Pt(),rapidity);
d7d7e825 1233 fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());
1234 fHistograms->FillHistogram("MC_Pi0_R", particle->R());
1235 fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());
1236 fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1237 fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
63e16c52 1238 if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_Fiducial", particle->Pt());
3c45d101 1239
1240 switch(GetProcessType(fGCMCEvent)){
1241 case kProcSD:
1242 fHistograms->FillHistogram("MC_SD_EvtQ3_Pi0_Pt", particle->Pt());
1243 break;
1244 case kProcDD:
1245 fHistograms->FillHistogram("MC_DD_EvtQ3_Pi0_Pt", particle->Pt());
1246 break;
1247 case kProcND:
1248 fHistograms->FillHistogram("MC_ND_EvtQ3_Pi0_Pt", particle->Pt());
1249 break;
1250 default:
1251 AliError("Unknown Process");
1252 }
63e16c52 1253
01b7fdcc 1254 if(gammaEtaCut && gammaRCut){
d7d7e825 1255 // if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
1256 fHistograms->FillHistogram("MC_Pi0_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1257 fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
63e16c52 1258 if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_withinAcceptance_Fiducial", particle->Pt());
1259
d7d7e825 1260 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1261 fHistograms->FillHistogram("MC_Pi0_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1262 fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1263 fHistograms->FillHistogram("MC_Pi0_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
10e3319b 1264 fHistograms->FillHistogram("MC_Pi0_ConvGamma_OpeningAngle_Pt", particle->Pt(),GetMCOpeningAngle(daughter0,daughter1));
1265 fHistograms->FillHistogram("MC_Pi0_ConvGamma_PtGamma_Pt", particle->Pt(),daughter0->Pt());
1266 fHistograms->FillHistogram("MC_Pi0_ConvGamma_PtGamma_Pt", particle->Pt(),daughter1->Pt());
1267
48682642 1268 Double_t alfa=0.;
77ac6f3e 1269 if((daughter0->Energy()+daughter1->Energy()) > 0.){
48682642 1270 alfa= TMath::Abs((daughter0->Energy()-daughter1->Energy())/(daughter0->Energy()+daughter1->Energy()));
1271 }
1272 fHistograms->FillHistogram("MC_Pi0_alpha",alfa);
1273 if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_ConvGamma_withinAcceptance_Fiducial", particle->Pt());
1274
d7d7e825 1275 }
1276 }
1277 }
1278 }
1279
1280 if(particle->GetPdgCode()==221){ //Eta
1281 fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());
1282 fHistograms->FillHistogram("MC_Eta_Rapid", rapidity);
1283 fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);
1284 fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());
dc2883e4 1285 fHistograms->FillHistogram("MC_Eta_Pt_vs_Rapid", particle->Pt(),rapidity);
d7d7e825 1286 fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());
1287 fHistograms->FillHistogram("MC_Eta_R", particle->R());
1288 fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());
1289 fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1290 fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
1291
01b7fdcc 1292 if(gammaEtaCut && gammaRCut){
a0b94e5c 1293 // if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
d7d7e825 1294 fHistograms->FillHistogram("MC_Eta_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1295 fHistograms->FillHistogram("MC_Eta_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
1296 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1297 fHistograms->FillHistogram("MC_Eta_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1298 fHistograms->FillHistogram("MC_Eta_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1299 fHistograms->FillHistogram("MC_Eta_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
10e3319b 1300 fHistograms->FillHistogram("MC_Eta_ConvGamma_OpeningAngle_Pt", particle->Pt(),GetMCOpeningAngle(daughter0,daughter1));
1301 fHistograms->FillHistogram("MC_Eta_ConvGamma_PtGamma_Pt", particle->Pt(),daughter0->Pt());
1302 fHistograms->FillHistogram("MC_Eta_ConvGamma_PtGamma_Pt", particle->Pt(),daughter1->Pt());
1303
d7d7e825 1304 }
1305
1306 }
1307
1308 }
1309
1310 // all motherparticles with 2 gammas as daughters
1311 fHistograms->FillHistogram("MC_Mother_R", particle->R());
1312 fHistograms->FillHistogram("MC_Mother_ZR", particle->Vz(),particle->R());
1313 fHistograms->FillHistogram("MC_Mother_XY", particle->Vx(),particle->Vy());
1314 fHistograms->FillHistogram("MC_Mother_Mass", particle->GetCalcMass());
1315 fHistograms->FillHistogram("MC_Mother_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1316 fHistograms->FillHistogram("MC_Mother_Energy", particle->Energy());
1317 fHistograms->FillHistogram("MC_Mother_Pt", particle->Pt());
1318 fHistograms->FillHistogram("MC_Mother_Eta", particle->Eta());
1319 fHistograms->FillHistogram("MC_Mother_Rapid", rapidity);
1320 fHistograms->FillHistogram("MC_Mother_Phi",tmpPhi);
1321 fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt",particle->GetMass(),particle->Pt());
a0b94e5c 1322
01b7fdcc 1323 if(gammaEtaCut && gammaRCut){
a0b94e5c 1324 // if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
d7d7e825 1325 fHistograms->FillHistogram("MC_Mother_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1326 fHistograms->FillHistogram("MC_Mother_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
1327 fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_withinAcceptance",particle->GetMass(),particle->Pt());
1328 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1329 fHistograms->FillHistogram("MC_Mother_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1330 fHistograms->FillHistogram("MC_Mother_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1331 fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_ConvGamma_withinAcceptance",particle->GetMass(),particle->Pt());
a0b94e5c 1332
d7d7e825 1333 }
1334
1335
1336 } // end passed R and eta cut
a0b94e5c 1337
d7d7e825 1338 } // end if(particle->GetNDaughters() == 2)
1339
1340 }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)
a0b94e5c 1341
d7d7e825 1342} // end ProcessMCData
1343
1344
1345
1346void AliAnalysisTaskGammaConversion::FillNtuple(){
1347 //Fills the ntuple with the different values
a0b94e5c 1348
d7d7e825 1349 if(fGammaNtuple == NULL){
1350 return;
1351 }
1352 Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
1353 for(Int_t i=0;i<numberOfV0s;i++){
1354 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};
1355 AliESDv0 * cV0 = fV0Reader->GetV0(i);
1356 Double_t negPID=0;
1357 Double_t posPID=0;
1358 fV0Reader->GetPIDProbability(negPID,posPID);
1359 values[0]=cV0->GetOnFlyStatus();
1360 values[1]=fV0Reader->CheckForPrimaryVertex();
1361 values[2]=negPID;
1362 values[3]=posPID;
1363 values[4]=fV0Reader->GetX();
1364 values[5]=fV0Reader->GetY();
1365 values[6]=fV0Reader->GetZ();
1366 values[7]=fV0Reader->GetXYRadius();
1367 values[8]=fV0Reader->GetMotherCandidateNDF();
1368 values[9]=fV0Reader->GetMotherCandidateChi2();
1369 values[10]=fV0Reader->GetMotherCandidateEnergy();
1370 values[11]=fV0Reader->GetMotherCandidateEta();
1371 values[12]=fV0Reader->GetMotherCandidatePt();
1372 values[13]=fV0Reader->GetMotherCandidateMass();
1373 values[14]=fV0Reader->GetMotherCandidateWidth();
1374 // values[15]=fV0Reader->GetMotherMCParticle()->Pt(); MOVED TO THE END, HAS TO BE CALLED AFTER HasSameMother NB: still has the same entry in the array
1375 values[16]=fV0Reader->GetOpeningAngle();
1376 values[17]=fV0Reader->GetNegativeTrackEnergy();
1377 values[18]=fV0Reader->GetNegativeTrackPt();
1378 values[19]=fV0Reader->GetNegativeTrackEta();
1379 values[20]=fV0Reader->GetNegativeTrackPhi();
1380 values[21]=fV0Reader->GetPositiveTrackEnergy();
1381 values[22]=fV0Reader->GetPositiveTrackPt();
1382 values[23]=fV0Reader->GetPositiveTrackEta();
1383 values[24]=fV0Reader->GetPositiveTrackPhi();
1384 values[25]=fV0Reader->HasSameMCMother();
1385 if(values[25] != 0){
1386 values[26]=fV0Reader->GetMotherMCParticlePDGCode();
1387 values[15]=fV0Reader->GetMotherMCParticle()->Pt();
1388 }
1389 fTotalNumberOfAddedNtupleEntries++;
1390 fGammaNtuple->Fill(values);
1391 }
1392 fV0Reader->ResetV0IndexNumber();
1393
1394}
1395
1396void AliAnalysisTaskGammaConversion::ProcessV0sNoCut(){
1397 // Process all the V0's without applying any cuts to it
a0b94e5c 1398
d7d7e825 1399 Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
1400 for(Int_t i=0;i<numberOfV0s;i++){
1401 /*AliESDv0 * cV0 = */fV0Reader->GetV0(i);
a0b94e5c 1402
d7d7e825 1403 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
cb90a330 1404 continue;
d7d7e825 1405 }
9640a3d1 1406
77880bd8 1407 // if( !fV0Reader->GetV0(i)->GetOnFlyStatus()){
1408 if( !fV0Reader->CheckV0FinderStatus(i)){
cb90a330 1409 continue;
9640a3d1 1410 }
1411
cb90a330 1412
1413 if( !((fV0Reader->GetNegativeESDTrack())->GetStatus() & AliESDtrack::kTPCrefit) ||
1414 !((fV0Reader->GetPositiveESDTrack())->GetStatus() & AliESDtrack::kTPCrefit) ){
1415 continue;
9640a3d1 1416 }
cb90a330 1417
ebcfaa7e 1418 if( fV0Reader->GetNegativeESDTrack()->GetSign()== fV0Reader->GetPositiveESDTrack()->GetSign()){
1419 continue;
1420 }
1421
1422 if( fV0Reader->GetNegativeESDTrack()->GetKinkIndex(0) > 0 ||
1423 fV0Reader->GetPositiveESDTrack()->GetKinkIndex(0) > 0) {
1424 continue;
1425 }
9c1cb6f7 1426 if(TMath::Abs(fV0Reader->GetMotherCandidateEta())> fV0Reader->GetEtaCut()){
1427 continue;
1428 }
1429 if(TMath::Abs(fV0Reader->GetPositiveTrackEta())> fV0Reader->GetEtaCut()){
1430 continue;
1431 }
1432 if(TMath::Abs(fV0Reader->GetNegativeTrackEta())> fV0Reader->GetEtaCut()){
1433 continue;
1434 }
1435 if((TMath::Abs(fV0Reader->GetZ())*fV0Reader->GetLineCutZRSlope())-fV0Reader->GetLineCutZValue() > fV0Reader->GetXYRadius() ){ // cuts out regions where we do not reconstruct
1436 continue;
1437 }
d7d7e825 1438 if(fDoMCTruth){
a0b94e5c 1439
d7d7e825 1440 if(fV0Reader->HasSameMCMother() == kFALSE){
1441 continue;
1442 }
a0b94e5c 1443
d7d7e825 1444 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1445 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
a0b94e5c 1446
d7d7e825 1447 if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1448 continue;
1449 }
26923b22 1450 if(negativeMC->GetPdgCode() == positiveMC->GetPdgCode()){
d7d7e825 1451 continue;
1452 }
a68437fb 1453
26923b22 1454 if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){ // id 5 is conversion
a68437fb 1455 continue;
1456 }
a0b94e5c 1457
d7d7e825 1458 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
a0b94e5c 1459
d7d7e825 1460 fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1461 fHistograms->FillHistogram("ESD_NoCutConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1462 fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1463 fHistograms->FillHistogram("ESD_NoCutConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1464 fHistograms->FillHistogram("ESD_NoCutConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1465 fHistograms->FillHistogram("ESD_NoCutConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1466 fHistograms->FillHistogram("ESD_NoCutConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1467 fHistograms->FillHistogram("ESD_NoCutConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1468 fHistograms->FillHistogram("ESD_NoCutConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1469 fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
a0b94e5c 1470
d7d7e825 1471 fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1472 fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
1473
1474 fHistograms->FillHistogram("ESD_NoCutConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1475 fHistograms->FillHistogram("ESD_NoCutConversion_R", fV0Reader->GetXYRadius());
1476 fHistograms->FillHistogram("ESD_NoCutConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1477 fHistograms->FillHistogram("ESD_NoCutConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
9640a3d1 1478 fHistograms->FillHistogram("ESD_NoCutConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1479 fHistograms->FillHistogram("ESD_NoCutConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1480 fHistograms->FillHistogram("ESD_NoCutConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1481 fHistograms->FillHistogram("ESD_NoCutConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1482
1483 fHistograms->FillHistogram("ESD_NoCutConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1484 fHistograms->FillHistogram("ESD_NoCutConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1485 fHistograms->FillHistogram("ESD_NoCutConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1486 fHistograms->FillHistogram("ESD_NoCutConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1487
d7d7e825 1488 //store MCTruth properties
1489 fHistograms->FillHistogram("ESD_NoCutConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
1490 fHistograms->FillHistogram("ESD_NoCutConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
1491 fHistograms->FillHistogram("ESD_NoCutConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
d7d7e825 1492 }
1493 }
1494 }
1495 fV0Reader->ResetV0IndexNumber();
1496}
1497
1498void AliAnalysisTaskGammaConversion::ProcessV0s(){
1499 // see header file for documentation
037dc2db 1500
1501
d7d7e825 1502 if(fWriteNtuple == kTRUE){
1503 FillNtuple();
1504 }
1505
1506 Int_t nSurvivingV0s=0;
5ce758b0 1507 fV0Reader->ResetNGoodV0s();
d7d7e825 1508 while(fV0Reader->NextV0()){
1509 nSurvivingV0s++;
1510
e60f3265 1511
1512 TVector3 vtxConv(fV0Reader->GetX(),fV0Reader->GetY(), fV0Reader->GetZ());
1513
d7d7e825 1514 //-------------------------- filling v0 information -------------------------------------
1515 fHistograms->FillHistogram("ESD_Conversion_R", fV0Reader->GetXYRadius());
1516 fHistograms->FillHistogram("ESD_Conversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1517 fHistograms->FillHistogram("ESD_Conversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1518 fHistograms->FillHistogram("ESD_Conversion_OpeningAngle", fV0Reader->GetOpeningAngle());
e60f3265 1519
1520 // Specific histograms for beam pipe studies
1521 if( TMath::Abs(fV0Reader->GetZ()) < fV0Reader->GetLineCutZValue() ){
1522 fHistograms->FillHistogram("ESD_Conversion_XY_BeamPipe", fV0Reader->GetX(),fV0Reader->GetY());
1523 fHistograms->FillHistogram("ESD_Conversion_RPhi_BeamPipe", vtxConv.Phi(),fV0Reader->GetXYRadius());
1524 }
1525
1526
d7d7e825 1527 fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());
1528 fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());
1529 fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());
1530 fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());
ebcfaa7e 1531 fHistograms->FillHistogram("ESD_E_nTPCClusters", fV0Reader->GetNegativeTracknTPCClusters());
1532 fHistograms->FillHistogram("ESD_E_nITSClusters", fV0Reader->GetNegativeTracknITSClusters());
5ce758b0 1533 if(fV0Reader->GetNegativeTracknTPCFClusters()!=0 && fV0Reader->GetNegativeTracknTPCClusters()!=0 ){
77ac6f3e 1534 Double_t eClsToF= (Double_t)fV0Reader->GetNegativeTracknTPCClusters()/(Double_t)fV0Reader->GetNegativeTracknTPCFClusters();
1535 fHistograms->FillHistogram("ESD_E_nTPCClustersToFP", fV0Reader->GetNegativeTrackP(),eClsToF );
5ce758b0 1536 fHistograms->FillHistogram("ESD_E_TPCchi2", fV0Reader->GetNegativeTrackTPCchi2()/(Double_t)fV0Reader->GetNegativeTracknTPCClusters());
1537 }
1538
1539
1540
d7d7e825 1541 fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());
1542 fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());
1543 fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());
1544 fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());
ebcfaa7e 1545 fHistograms->FillHistogram("ESD_P_nTPCClusters", fV0Reader->GetPositiveTracknTPCClusters());
1546 fHistograms->FillHistogram("ESD_P_nITSClusters", fV0Reader->GetPositiveTracknITSClusters());
5ce758b0 1547 if(fV0Reader->GetPositiveTracknTPCFClusters()!=0 && (Double_t)fV0Reader->GetPositiveTracknTPCClusters()!=0 ){
77ac6f3e 1548 Double_t pClsToF= (Double_t)fV0Reader->GetPositiveTracknTPCClusters()/(Double_t)fV0Reader->GetPositiveTracknTPCFClusters();
1549 fHistograms->FillHistogram("ESD_P_nTPCClustersToFP",fV0Reader->GetPositiveTrackP(), pClsToF);
5ce758b0 1550 fHistograms->FillHistogram("ESD_P_TPCchi2", fV0Reader->GetPositiveTrackTPCchi2()/(Double_t)fV0Reader->GetPositiveTracknTPCClusters());
1551 }
1552
1553
d7d7e825 1554 fHistograms->FillHistogram("ESD_ConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1555 fHistograms->FillHistogram("ESD_ConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1556 fHistograms->FillHistogram("ESD_ConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1557 fHistograms->FillHistogram("ESD_ConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1558 fHistograms->FillHistogram("ESD_ConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1559 fHistograms->FillHistogram("ESD_ConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1560 fHistograms->FillHistogram("ESD_ConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1561 fHistograms->FillHistogram("ESD_ConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1562 fHistograms->FillHistogram("ESD_ConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1563 fHistograms->FillHistogram("ESD_ConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1564
1565 fHistograms->FillHistogram("ESD_ConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1566 fHistograms->FillHistogram("ESD_ConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
9640a3d1 1567
1568 fHistograms->FillHistogram("ESD_ConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1569 fHistograms->FillHistogram("ESD_ConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1570 fHistograms->FillHistogram("ESD_ConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1571 fHistograms->FillHistogram("ESD_ConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1572
1573 fHistograms->FillHistogram("ESD_ConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1574 fHistograms->FillHistogram("ESD_ConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1575 fHistograms->FillHistogram("ESD_ConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1576 fHistograms->FillHistogram("ESD_ConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1577
9c1cb6f7 1578 Double_t negPID=0;
1579 Double_t posPID=0;
1580 fV0Reader->GetPIDProbability(negPID,posPID);
1581 fHistograms->FillHistogram("ESD_ConvGamma_E_EProbP",fV0Reader->GetNegativeTrackP(),negPID);
1582 fHistograms->FillHistogram("ESD_ConvGamma_P_EProbP",fV0Reader->GetPositiveTrackP(),posPID);
1583
1584 Double_t negPIDmupi=0;
1585 Double_t posPIDmupi=0;
1586 fV0Reader->GetPIDProbabilityMuonPion(negPIDmupi,posPIDmupi);
1587 fHistograms->FillHistogram("ESD_ConvGamma_E_mupiProbP",fV0Reader->GetNegativeTrackP(),negPIDmupi);
1588 fHistograms->FillHistogram("ESD_ConvGamma_P_mupiProbP",fV0Reader->GetPositiveTrackP(),posPIDmupi);
1589
48682642 1590 Double_t armenterosQtAlfa[2];
1591 fV0Reader->GetArmenterosQtAlfa(fV0Reader-> GetNegativeKFParticle(),
1592 fV0Reader-> GetPositiveKFParticle(),
1593 fV0Reader->GetMotherCandidateKFCombination(),
1594 armenterosQtAlfa);
1595
1596 fHistograms->FillHistogram("ESD_ConvGamma_alfa_qt",armenterosQtAlfa[1],armenterosQtAlfa[0]);
1597
037dc2db 1598
d7d7e825 1599 // begin mapping
1600 Int_t rBin = fHistograms->GetRBin(fV0Reader->GetXYRadius());
9640a3d1 1601 Int_t zBin = fHistograms->GetZBin(fV0Reader->GetZ());
d7d7e825 1602 Int_t phiBin = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());
9c1cb6f7 1603 Double_t rFMD=30;
64f4118c 1604 Double_t rITSTPCMin=50;
1605 Double_t rITSTPCMax=80;
e60f3265 1606
9640a3d1 1607
ebcfaa7e 1608 // Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();
d7d7e825 1609
1610 TString nameESDMappingPhiR="";
e158cbc3 1611 nameESDMappingPhiR.Form("ESD_Conversion_Mapping_Phi%02d_R%02d",phiBin,rBin);
ebcfaa7e 1612 //fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);
d7d7e825 1613
1614 TString nameESDMappingPhi="";
e158cbc3 1615 nameESDMappingPhi.Form("ESD_Conversion_Mapping_Phi%02d",phiBin);
ebcfaa7e 1616 //fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);
d7d7e825 1617
1618 TString nameESDMappingR="";
e158cbc3 1619 nameESDMappingR.Form("ESD_Conversion_Mapping_R%02d",rBin);
ebcfaa7e 1620 //fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);
d7d7e825 1621
1622 TString nameESDMappingPhiInR="";
e158cbc3 1623 nameESDMappingPhiInR.Form("ESD_Conversion_Mapping_Phi_in_R_%02d",rBin);
9640a3d1 1624 // fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
1625 fHistograms->FillHistogram(nameESDMappingPhiInR, vtxConv.Phi());
1626
1627 TString nameESDMappingZInR="";
1628 nameESDMappingZInR.Form("ESD_Conversion_Mapping_Z_in_R_%02d",rBin);
1629 fHistograms->FillHistogram(nameESDMappingZInR, fV0Reader->GetZ());
1630
1631 TString nameESDMappingPhiInZ="";
1632 nameESDMappingPhiInZ.Form("ESD_Conversion_Mapping_Phi_in_Z_%02d",zBin);
1633 // fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
1634 fHistograms->FillHistogram(nameESDMappingPhiInZ, vtxConv.Phi());
1635
9c1cb6f7 1636 if(fV0Reader->GetXYRadius()<rFMD){
1637 TString nameESDMappingFMDPhiInZ="";
1638 nameESDMappingFMDPhiInZ.Form("ESD_Conversion_Mapping_FMD_Phi_in_Z_%02d",zBin);
1639 fHistograms->FillHistogram(nameESDMappingFMDPhiInZ, vtxConv.Phi());
1640 }
1641
64f4118c 1642 if(fV0Reader->GetXYRadius()>rITSTPCMin && fV0Reader->GetXYRadius()<rITSTPCMax){
1643 TString nameESDMappingITSTPCPhiInZ="";
1644 nameESDMappingITSTPCPhiInZ.Form("ESD_Conversion_Mapping_ITSTPC_Phi_in_Z_%02d",zBin);
1645 fHistograms->FillHistogram(nameESDMappingITSTPCPhiInZ, vtxConv.Phi());
1646 }
9c1cb6f7 1647
9640a3d1 1648 TString nameESDMappingRInZ="";
1649 nameESDMappingRInZ.Form("ESD_Conversion_Mapping_R_in_Z_%02d",zBin);
1650 fHistograms->FillHistogram(nameESDMappingRInZ, fV0Reader->GetXYRadius());
1651
1652 if(fV0Reader->GetMotherCandidatePt() > fLowPtMapping && fV0Reader->GetMotherCandidatePt()< fHighPtMapping){
1653 TString nameESDMappingMidPtPhiInR="";
1654 nameESDMappingMidPtPhiInR.Form("ESD_Conversion_Mapping_MidPt_Phi_in_R_%02d",rBin);
1655 fHistograms->FillHistogram(nameESDMappingMidPtPhiInR, vtxConv.Phi());
1656
1657 TString nameESDMappingMidPtZInR="";
1658 nameESDMappingMidPtZInR.Form("ESD_Conversion_Mapping_MidPt_Z_in_R_%02d",rBin);
1659 fHistograms->FillHistogram(nameESDMappingMidPtZInR, fV0Reader->GetZ());
1660
1661 TString nameESDMappingMidPtPhiInZ="";
1662 nameESDMappingMidPtPhiInZ.Form("ESD_Conversion_Mapping_MidPt_Phi_in_Z_%02d",zBin);
1663 fHistograms->FillHistogram(nameESDMappingMidPtPhiInZ, vtxConv.Phi());
9c1cb6f7 1664 if(fV0Reader->GetXYRadius()<rFMD){
1665 TString nameESDMappingMidPtFMDPhiInZ="";
1666 nameESDMappingMidPtFMDPhiInZ.Form("ESD_Conversion_Mapping_MidPt_FMD_Phi_in_Z_%02d",zBin);
1667 fHistograms->FillHistogram(nameESDMappingMidPtFMDPhiInZ, vtxConv.Phi());
1668 }
1669
9640a3d1 1670
1671 TString nameESDMappingMidPtRInZ="";
1672 nameESDMappingMidPtRInZ.Form("ESD_Conversion_Mapping_MidPt_R_in_Z_%02d",zBin);
1673 fHistograms->FillHistogram(nameESDMappingMidPtRInZ, fV0Reader->GetXYRadius());
1674
1675 }
1676
1677
d7d7e825 1678 // end mapping
1679
6c84d371 1680 new((*fKFReconstructedGammasTClone)[fKFReconstructedGammasTClone->GetEntriesFast()]) AliKFParticle(*fV0Reader->GetMotherCandidateKFCombination());
037dc2db 1681 fKFReconstructedGammasV0Index.push_back(fV0Reader->GetCurrentV0IndexNumber()-1);
6c84d371 1682 // fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());
d7d7e825 1683 fElectronv1.push_back(fV0Reader->GetCurrentV0()->GetPindex());
1684 fElectronv2.push_back(fV0Reader->GetCurrentV0()->GetNindex());
a0b94e5c 1685
d7d7e825 1686 //----------------------------------- checking for "real" conversions (MC match) --------------------------------------
1687 if(fDoMCTruth){
64f4118c 1688 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1689 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
d7d7e825 1690
a68437fb 1691 if(fV0Reader->HasSameMCMother() == kFALSE){
64f4118c 1692 fHistograms->FillHistogram("ESD_TrueConvCombinatorial_R", fV0Reader->GetXYRadius());
1693 if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
1694 fHistograms->FillHistogram("ESD_TrueConvCombinatorialElec_R", fV0Reader->GetXYRadius());
1695 }
a68437fb 1696 continue;
1697 }
64f4118c 1698 // Moved up to check true electron background
1699 // TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1700 // TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
a68437fb 1701
1702 if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1703 continue;
1704 }
1705
1706 if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
1707 continue;
1708 }
8de8695b 1709
1710 UInt_t statusPos = fV0Reader->GetPositiveESDTrack()->GetStatus();
1711 UInt_t statusNeg = fV0Reader->GetNegativeESDTrack()->GetStatus();
1712 UChar_t itsPixelPos = fV0Reader->GetPositiveESDTrack()->GetITSClusterMap();
1713 UChar_t itsPixelNeg = fV0Reader->GetNegativeESDTrack()->GetITSClusterMap();
1714
1715 // Using the UniqueID Phojet does not get the Dalitz right
1716 // if( (negativeMC->GetUniqueID() == 4 && positiveMC->GetUniqueID() ==4) ||
1717 // (negativeMC->GetUniqueID() == 0 && positiveMC->GetUniqueID() ==0) ){// fill r distribution for Dalitz decays
1718 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 111){ //pi0
1719 fHistograms->FillHistogram("ESD_TrueDalitzContamination_R", fV0Reader->GetXYRadius());
1720 fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_R", fV0Reader->GetXYRadius());
1721 //--------Histos for HFE
1722
1723 if(statusPos & AliESDtrack::kTOFpid){
1724 fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_SinglePos_R", fV0Reader->GetXYRadius());
1725 if( TESTBIT(itsPixelPos, 0) ){
1726 fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_SinglePos_kFirst_R", fV0Reader->GetXYRadius());
1727 }
037dc2db 1728 }
8de8695b 1729 if(statusNeg & AliESDtrack::kTOFpid){
1730 fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_SingleNeg_R", fV0Reader->GetXYRadius());
1731 if( TESTBIT(itsPixelNeg, 0) ){
1732 fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_SingleNeg_kFirst_R", fV0Reader->GetXYRadius());
1733 }
64f4118c 1734 }
8de8695b 1735 //--------------------------------------------------------
64f4118c 1736
037dc2db 1737 }
8de8695b 1738 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 221){ //eta
1739 fHistograms->FillHistogram("ESD_TrueConvDalitzEta_R", fV0Reader->GetXYRadius());
1740 }
1741
1742 //}
a68437fb 1743
1744 if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){// check if the daughters come from a conversion
1745 continue;
1746 }
a0b94e5c 1747
d7d7e825 1748 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
a68437fb 1749
26923b22 1750 if(fDoCF){
1751 Double_t containerInput[3];
1752 containerInput[0] = fV0Reader->GetMotherCandidatePt();
1753 containerInput[1] = fV0Reader->GetMotherCandidateEta();
1754 containerInput[2] = fV0Reader->GetMotherCandidateMass();
1755 fCFManager->GetParticleContainer()->Fill(containerInput,kStepTrueGamma); // for CF
1756 }
1757
d7d7e825 1758 fHistograms->FillHistogram("ESD_TrueConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1759 fHistograms->FillHistogram("ESD_TrueConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1760 fHistograms->FillHistogram("ESD_TrueConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1761 fHistograms->FillHistogram("ESD_TrueConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1762 fHistograms->FillHistogram("ESD_TrueConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1763 fHistograms->FillHistogram("ESD_TrueConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1764 fHistograms->FillHistogram("ESD_TrueConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1765 fHistograms->FillHistogram("ESD_TrueConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1766 fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1767 fHistograms->FillHistogram("ESD_TrueConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1768 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters());
1769 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters());
1770 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters(),fV0Reader->GetMotherCandidateMass());
1771 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters(),fV0Reader->GetMotherCandidateMass());
a0b94e5c 1772
d7d7e825 1773 fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1774 fHistograms->FillHistogram("ESD_TrueConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
a0b94e5c 1775
d7d7e825 1776
1777 fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1778 fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());
1779 fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1780 fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
9640a3d1 1781
8de8695b 1782 //----Histos for HFE--------------------------------------
1783
1784 if(statusPos & AliESDtrack::kTOFpid){
1785 fHistograms->FillHistogram("ESD_TrueConversion_SinglePos_R", positiveMC->R(),fV0Reader->GetPositiveMCParticle()->Pt());
1786 if( TESTBIT(itsPixelPos, 0) ){
1787 fHistograms->FillHistogram("ESD_TrueConversion_SinglePos_kFirst_R", positiveMC->R(),fV0Reader->GetPositiveMCParticle()->Pt());
1788 }
1789 }
1790 if(statusNeg & AliESDtrack::kTOFpid){
1791 fHistograms->FillHistogram("ESD_TrueConversion_SingleNeg_R", negativeMC->R(),fV0Reader->GetNegativeMCParticle()->Pt());
1792 if( TESTBIT(itsPixelNeg, 0) ){
1793 fHistograms->FillHistogram("ESD_TrueConversion_SingleNeg_kFirst_R", negativeMC->R(),fV0Reader->GetNegativeMCParticle()->Pt());
1794 }
1795 }
1796 //--------------------------------------------------------
1797
9640a3d1 1798 fHistograms->FillHistogram("ESD_TrueConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1799 fHistograms->FillHistogram("ESD_TrueConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1800 fHistograms->FillHistogram("ESD_TrueConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1801 fHistograms->FillHistogram("ESD_TrueConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
96ade8ef 1802 if (fV0Reader->GetMotherCandidateP() != 0) {
1803 fHistograms->FillHistogram("ESD_TrueConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1804 fHistograms->FillHistogram("ESD_TrueConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1805 } else { cout << "Error::fV0Reader->GetNegativeTrackP() == 0 !!!" << endl; }
1806
9640a3d1 1807 fHistograms->FillHistogram("ESD_TrueConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1808 fHistograms->FillHistogram("ESD_TrueConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
70ef88b5 1809 fHistograms->FillHistogram("ESD_TrueConvGamma_alfa_qt",armenterosQtAlfa[1],armenterosQtAlfa[0]);
1810
9640a3d1 1811
a0b94e5c 1812
d7d7e825 1813 //store MCTruth properties
1814 fHistograms->FillHistogram("ESD_TrueConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
1815 fHistograms->FillHistogram("ESD_TrueConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
1816 fHistograms->FillHistogram("ESD_TrueConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
48682642 1817
d7d7e825 1818 //resolution
1819 Double_t mcpt = fV0Reader->GetMotherMCParticle()->Pt();
1820 Double_t esdpt = fV0Reader->GetMotherCandidatePt();
9c1cb6f7 1821 Double_t resdPt = 0.;
4a6157dc 1822 if(mcpt > 0){
9c1cb6f7 1823 resdPt = ((esdpt - mcpt)/mcpt)*100.;
d7d7e825 1824 }
4a6157dc 1825 else if(mcpt < 0){
1826 cout<<"Pt of MC particle is negative, this will cause wrong calculation of resPt"<<endl;
1827 }
d7d7e825 1828
ca6d4600 1829 fHistograms->FillHistogram("Resolution_Gamma_dPt_Pt", mcpt, resdPt);
d7d7e825 1830 fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);
1831 fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);
ca6d4600 1832 fHistograms->FillHistogram("Resolution_Gamma_dPt_Phi", fV0Reader->GetMotherCandidatePhi(), resdPt);
d7d7e825 1833
9c1cb6f7 1834 Double_t resdZ = 0.;
d7d7e825 1835 if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){
9c1cb6f7 1836 resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100.;
d7d7e825 1837 }
9c1cb6f7 1838 Double_t resdZAbs = 0.;
037dc2db 1839 resdZAbs = (fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz());
1840
1841 fHistograms->FillHistogram("Resolution_dZAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdZAbs);
d7d7e825 1842 fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);
1843 fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());
1844 fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());
ca6d4600 1845
1846 // new for dPt_Pt-histograms for Electron and Positron
1847 Double_t mcEpt = fV0Reader->GetNegativeMCParticle()->Pt();
9c1cb6f7 1848 Double_t resEdPt = 0.;
ca6d4600 1849 if (mcEpt > 0){
9c1cb6f7 1850 resEdPt = ((fV0Reader->GetNegativeTrackPt()-mcEpt)/mcEpt)*100.;
ca6d4600 1851 }
1852 UInt_t statusN = fV0Reader->GetNegativeESDTrack()->GetStatus();
1853 // AliESDtrack * negTrk = fV0Reader->GetNegativeESDTrack();
1854 UInt_t kTRDoutN = (statusN & AliESDtrack::kTRDout);
1855
77ac6f3e 1856 Int_t nITSclsE= fV0Reader->GetNegativeTracknITSClusters();
ca6d4600 1857 // filling Resolution_Pt_dPt with respect to the Number of ITS clusters for Positrons
77ac6f3e 1858 switch(nITSclsE){
ca6d4600 1859 case 0: // 0 ITS clusters
1860 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS0", mcEpt, resEdPt);
1861 break;
1862 case 1: // 1 ITS cluster
1863 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS1", mcEpt, resEdPt);
1864 break;
1865 case 2: // 2 ITS clusters
1866 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS2", mcEpt, resEdPt);
1867 break;
1868 case 3: // 3 ITS clusters
1869 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS3", mcEpt, resEdPt);
1870 break;
1871 case 4: // 4 ITS clusters
1872 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS4", mcEpt, resEdPt);
1873 break;
1874 case 5: // 5 ITS clusters
1875 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS5", mcEpt, resEdPt);
1876 break;
1877 case 6: // 6 ITS clusters
1878 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS6", mcEpt, resEdPt);
1879 break;
1880 }
1881 //Filling histograms with respect to Electron resolution
1882 fHistograms->FillHistogram("Resolution_E_dPt_Pt", mcEpt, resEdPt);
1883 fHistograms->FillHistogram("Resolution_E_dPt_Phi", fV0Reader->GetNegativeTrackPhi(), resEdPt);
1884 if(kTRDoutN){
1885 fHistograms->FillHistogram("Resolution_E_nTRDtracklets_ESDPt", fV0Reader->GetNegativeTrackPt(), fV0Reader->GetNegativeESDTrack()->GetTRDntracklets());
1886 fHistograms->FillHistogram("Resolution_E_nTRDtracklets_MCPt", mcEpt, fV0Reader->GetNegativeESDTrack()->GetTRDntracklets());
1887 fHistograms->FillHistogram("Resolution_E_nTRDclusters_ESDPt",fV0Reader->GetNegativeTrackPt(), fV0Reader->GetNegativeESDTrack()->GetTRDncls());
1888 fHistograms->FillHistogram("Resolution_E_nTRDclusters_MCPt",mcEpt, fV0Reader->GetNegativeESDTrack()->GetTRDncls());
1889 fHistograms->FillHistogram("Resolution_E_TRDsignal_ESDPt", fV0Reader->GetNegativeTrackPt(), fV0Reader->GetNegativeESDTrack()->GetTRDsignal());
1890 }
1891
1892 Double_t mcPpt = fV0Reader->GetPositiveMCParticle()->Pt();
1893 Double_t resPdPt = 0;
1894 if (mcPpt > 0){
9c1cb6f7 1895 resPdPt = ((fV0Reader->GetPositiveTrackPt()-mcPpt)/mcPpt)*100.;
ca6d4600 1896 }
1897
1898 UInt_t statusP = fV0Reader->GetPositiveESDTrack()->GetStatus();
1899// AliESDtrack * posTr= fV0Reader->GetPositiveESDTrack();
1900 UInt_t kTRDoutP = (statusP & AliESDtrack::kTRDout);
1901
77ac6f3e 1902 Int_t nITSclsP = fV0Reader->GetPositiveTracknITSClusters();
ca6d4600 1903 // filling Resolution_Pt_dPt with respect to the Number of ITS clusters for Positrons
77ac6f3e 1904 switch(nITSclsP){
ca6d4600 1905 case 0: // 0 ITS clusters
1906 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS0", mcPpt, resPdPt);
1907 break;
1908 case 1: // 1 ITS cluster
1909 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS1", mcPpt, resPdPt);
1910 break;
1911 case 2: // 2 ITS clusters
1912 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS2", mcPpt, resPdPt);
1913 break;
1914 case 3: // 3 ITS clusters
1915 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS3", mcPpt, resPdPt);
1916 break;
1917 case 4: // 4 ITS clusters
1918 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS4", mcPpt, resPdPt);
1919 break;
1920 case 5: // 5 ITS clusters
1921 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS5", mcPpt, resPdPt);
1922 break;
1923 case 6: // 6 ITS clusters
1924 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS6", mcPpt, resPdPt);
1925 break;
1926 }
1927 //Filling histograms with respect to Positron resolution
1928 fHistograms->FillHistogram("Resolution_P_dPt_Pt", mcPpt, resPdPt);
1929 fHistograms->FillHistogram("Resolution_P_dPt_Phi", fV0Reader->GetPositiveTrackPhi(), resPdPt);
1930 if(kTRDoutP){
1931 fHistograms->FillHistogram("Resolution_P_nTRDtracklets_ESDPt", fV0Reader->GetPositiveTrackPt(), fV0Reader->GetPositiveESDTrack()->GetTRDntracklets());
1932 fHistograms->FillHistogram("Resolution_P_nTRDtracklets_MCPt", mcPpt, fV0Reader->GetPositiveESDTrack()->GetTRDntracklets());
1933 fHistograms->FillHistogram("Resolution_P_nTRDclusters_ESDPt",fV0Reader->GetPositiveTrackPt(), fV0Reader->GetPositiveESDTrack()->GetTRDncls());
1934 fHistograms->FillHistogram("Resolution_P_nTRDclusters_MCPt",mcPpt, fV0Reader->GetPositiveESDTrack()->GetTRDncls());
1935 fHistograms->FillHistogram("Resolution_P_TRDsignal_ESDPt", fV0Reader->GetPositiveTrackPt(), fV0Reader->GetPositiveESDTrack()->GetTRDsignal());
1936 }
1937
1938
9c1cb6f7 1939 Double_t resdR = 0.;
d7d7e825 1940 if(fV0Reader->GetNegativeMCParticle()->R() != 0){
9c1cb6f7 1941 resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100.;
d7d7e825 1942 }
9c1cb6f7 1943 Double_t resdRAbs = 0.;
037dc2db 1944 resdRAbs = (fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R());
1945
1946 fHistograms->FillHistogram("Resolution_dRAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdRAbs);
d7d7e825 1947 fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);
1948 fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());
1949 fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());
ca6d4600 1950 fHistograms->FillHistogram("Resolution_R_dPt", fV0Reader->GetNegativeMCParticle()->R(), resdPt);
1951
9c1cb6f7 1952 Double_t resdPhiAbs=0.;
1953 resdPhiAbs=0.;
037dc2db 1954 resdPhiAbs= (fV0Reader->GetMotherCandidatePhi()-fV0Reader->GetNegativeMCParticle()->Phi());
1955 fHistograms->FillHistogram("Resolution_dPhiAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdPhiAbs);
1956
d7d7e825 1957 }//if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22)
d7d7e825 1958 }//if(fDoMCTruth)
1959 }//while(fV0Reader->NextV0)
1960 fHistograms->FillHistogram("ESD_NumberOfSurvivingV0s", nSurvivingV0s);
1961 fHistograms->FillHistogram("ESD_NumberOfV0s", fV0Reader->GetNumberOfV0s());
b5832f95 1962 fHistograms->FillHistogram("ESD_NumberOfContributorsVtx", fV0Reader->GetNumberOfContributorsVtx());
cb90a330 1963
1964 fV0Reader->ResetV0IndexNumber();
d7d7e825 1965}
1966
6441e967 1967void AliAnalysisTaskGammaConversion::AddToAODBranch(TClonesArray * branch, AliAODPWG4Particle & particle) {
1968 //See header file for documentation
1969
1970 Int_t i = branch->GetEntriesFast();
1971 if(! (fOutputAODClassName.Contains("Correlation")) ) {
1972 new((*branch)[i]) AliAODPWG4Particle(particle);
1973 } else {
1974 new((*branch)[i]) AliAODPWG4ParticleCorrelation(particle);
1975 }
1976}
1977
1978void AliAnalysisTaskGammaConversion::AddToAODBranch(TClonesArray * branch, AliGammaConversionAODObject & particle) {
1979 //See header file for documentation
1980
1981 Int_t i = branch->GetEntriesFast();
1982 new((*branch)[i]) AliGammaConversionAODObject(particle);
1983}
1984
1985void AliAnalysisTaskGammaConversion::AddToAODBranch(TClonesArray * branch, AliAODConversionParticle & particle) {
1986 //See header file for documentation
1987
1988 Int_t i = branch->GetEntriesFast();
1989 new((*branch)[i]) AliAODConversionParticle(particle);
1990}
1991
1992
d7d7e825 1993void AliAnalysisTaskGammaConversion::FillAODWithConversionGammas(){
1994 // Fill AOD with reconstructed Gamma
a0b94e5c 1995
6c84d371 1996 for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
332f1f44 1997 AliKFParticle * gammakf = dynamic_cast<AliKFParticle*>(fKFReconstructedGammasTClone->At(gammaIndex));
6441e967 1998
1999 if(fOutputAODClassName.Contains("AliAODPWG4Particle")) {
2000 AliAODPWG4Particle gamma = AliAODPWG4Particle(gammakf->Px(),gammakf->Py(),gammakf->Pz(), gammakf->E());
2001 //gamma.SetLabel(-1);//How to get the MC label of the reconstructed gamma?
2002 gamma.SetTrackLabel( fElectronv1[gammaIndex], fElectronv2[gammaIndex] ); //How to get the MC label of the 2 electrons that form the gamma?
2003 gamma.SetDetector("CTS"); //tag the gamma as reconstructed in the central barrel
2004 gamma.SetPdg(AliPID::kEleCon); //photon id
2005 gamma.SetTag(-1); //Here I usually put a flag saying that montecarlo says it is prompt, decay fragmentation photon, or hadrons or whatever
caf688bd 2006 //PH gamma.SetChi2(gammakf->Chi2());
6441e967 2007
2008 AddToAODBranch(fAODGamma, gamma);
2009
2010 } else if(fOutputAODClassName.Contains("ConversionParticle")) {
2011 TLorentzVector momentum(gammakf->Px(),gammakf->Py(),gammakf->Pz(), gammakf->E());
2012 AliAODConversionParticle gamma = AliAODConversionParticle(momentum);
2013 //gamma.SetLabel(-1);//How to get the MC label of the reconstructed gamma?
2014 gamma.SetTrackLabels( fElectronv1[gammaIndex], fElectronv2[gammaIndex] ); //How to get the MC label of the 2 electrons that form the gamma?
2015 //gamma.SetPdg(AliPID::kEleCon); //photon id
2016 //gamma.SetTag(-1); //Here I usually put a flag saying that montecarlo says it is prompt, decay fragmentation photon, or hadrons or whatever
caf688bd 2017 //PH gamma.SetChi2(gammakf->Chi2());
6441e967 2018 gamma.SetTrackLabels( fElectronv1[gammaIndex], fElectronv2[gammaIndex] );
2019 gamma.SetESDEvent(dynamic_cast<AliESDEvent*>(InputEvent()));
2020 AddToAODBranch(fAODGamma, gamma);
2021
2022
2023
2024 } else {
2025 AliGammaConversionAODObject gamma;
2026 gamma.SetPx(gammakf->GetPx());
2027 gamma.SetPy(gammakf->GetPy());
2028 gamma.SetPz(gammakf->GetPz());
2029 gamma.SetLabel1(fElectronv1[gammaIndex]);
2030 gamma.SetLabel2(fElectronv2[gammaIndex]);
2031 gamma.SetChi2(gammakf->Chi2());
2032 gamma.SetE(gammakf->E());
2033 gamma.SetESDEvent(dynamic_cast<AliESDEvent*>(InputEvent()));
2034 AddToAODBranch(fAODGamma, gamma);
2035 }
a0b94e5c 2036
6441e967 2037 }
d7d7e825 2038}
6272370b 2039void AliAnalysisTaskGammaConversion::ProcessGammasForOmegaMesonAnalysis(){
2040 // omega meson analysis pi0+gamma decay
2041 for(Int_t firstPi0Index=0;firstPi0Index<fKFReconstructedPi0sTClone->GetEntriesFast();firstPi0Index++){
48682642 2042 AliKFParticle * omegaCandidatePi0Daughter = (AliKFParticle *)fKFReconstructedPi0sTClone->At(firstPi0Index);
6272370b 2043 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
48682642 2044
6272370b 2045 AliKFParticle * omegaCandidateGammaDaughter = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2046 if(fGammav1[firstPi0Index]==firstGammaIndex || fGammav2[firstPi0Index]==firstGammaIndex){
2047 continue;
2048 }
2049
48682642 2050 AliKFParticle omegaCandidate(*omegaCandidatePi0Daughter,*omegaCandidateGammaDaughter);
6272370b 2051 Double_t massOmegaCandidate = 0.;
2052 Double_t widthOmegaCandidate = 0.;
2053
48682642 2054 omegaCandidate.GetMass(massOmegaCandidate,widthOmegaCandidate);
6272370b 2055
e9aea39f 2056 if ( massOmegaCandidate > 733 && massOmegaCandidate < 833 ) {
6441e967 2057 //AddOmegaToAOD(&omegaCandidate, massOmegaCandidate, firstPi0Index, firstGammaIndex);
c59360eb 2058 }
2059
48682642 2060 fHistograms->FillHistogram("ESD_Omega_InvMass_vs_Pt",massOmegaCandidate ,omegaCandidate.GetPt());
6272370b 2061 fHistograms->FillHistogram("ESD_Omega_InvMass",massOmegaCandidate);
48682642 2062
2063 //delete omegaCandidate;
2064
2065 }// end of omega reconstruction in pi0+gamma channel
2066
2067 if(fDoJet == kTRUE){
2068 AliKFParticle* negPiKF=NULL;
2069 AliKFParticle* posPiKF=NULL;
2070
2071 // look at the pi+pi+pi0 channel
2072 for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
2073 AliESDtrack* posTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
2074 if (posTrack->GetSign()<0) continue;
c0d9051e 2075 if(TMath::Abs(fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kPion))>2.) continue;
48682642 2076 if (posPiKF) delete posPiKF; posPiKF=NULL;
2077 posPiKF = new AliKFParticle( *(posTrack) ,211);
2078
2079 for(Int_t jCh=0;jCh<fChargedParticles->GetEntriesFast();jCh++){
2080 AliESDtrack* negTrack = (AliESDtrack*)(fChargedParticles->At(jCh));
2081 if( negTrack->GetSign()>0) continue;
c0d9051e 2082 if(TMath::Abs(fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kPion))>2.) continue;
48682642 2083 if (negPiKF) delete negPiKF; negPiKF=NULL;
2084 negPiKF = new AliKFParticle( *(negTrack) ,-211);
2085 AliKFParticle omegaCandidatePipPinPi0(*omegaCandidatePi0Daughter,*posPiKF,*negPiKF);
2086 Double_t massOmegaCandidatePipPinPi0 = 0.;
2087 Double_t widthOmegaCandidatePipPinPi0 = 0.;
2088
2089 omegaCandidatePipPinPi0.GetMass(massOmegaCandidatePipPinPi0,widthOmegaCandidatePipPinPi0);
ef2e2748 2090
2091 if ( massOmegaCandidatePipPinPi0 > 733 && massOmegaCandidatePipPinPi0 < 833 ) {
6441e967 2092 // AddOmegaToAOD(&omegaCandidatePipPinPi0, massOmegaCandidatePipPinPi0, -1, -1);
ef2e2748 2093 }
48682642 2094
2095 fHistograms->FillHistogram("ESD_OmegaPipPinPi0_InvMass_vs_Pt",massOmegaCandidatePipPinPi0 ,omegaCandidatePipPinPi0.GetPt());
2096 fHistograms->FillHistogram("ESD_OmegaPipPinPi0_InvMass",massOmegaCandidatePipPinPi0);
2097
2098 // delete omegaCandidatePipPinPi0;
2099 }
2100 }
48682642 2101
45936300 2102 if (posPiKF) delete posPiKF; posPiKF=NULL; if (negPiKF) delete negPiKF; negPiKF=NULL;
48682642 2103
45936300 2104 } // checking ig gammajet because in that case the chargedparticle list is created
6272370b 2105
6272370b 2106 }
48682642 2107
2108 if(fCalculateBackground){
5ce758b0 2109
2110 AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
2111
2112 Int_t zbin= bgHandler->GetZBinIndex(fV0Reader->GetVertexZ());
2113 Int_t mbin = 0;
2114 if(fUseTrackMultiplicityForBG == kTRUE){
2115 mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2116 }
2117 else{
2118 mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->GetNGoodV0s());
2119 }
2120
2121 AliGammaConversionBGHandler::GammaConversionVertex *bgEventVertex = NULL;
2122
48682642 2123 // Background calculation for the omega
2124 for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
5ce758b0 2125 AliGammaConversionKFVector * previousEventV0s = bgHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);
2126
2127 if(fMoveParticleAccordingToVertex == kTRUE){
2128 bgEventVertex = bgHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
2129 }
48682642 2130 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2131 AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
5ce758b0 2132
2133 if(fMoveParticleAccordingToVertex == kTRUE){
2134 MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
2135 }
2136
48682642 2137 for(Int_t firstPi0Index=0;firstPi0Index<fKFReconstructedPi0sTClone->GetEntriesFast();firstPi0Index++){
2138 AliKFParticle * omegaCandidatePi0Daughter = (AliKFParticle *)fKFReconstructedPi0sTClone->At(firstPi0Index);
2139 AliKFParticle * omegaBckCandidate = new AliKFParticle(*omegaCandidatePi0Daughter,previousGoodV0);
2140 Double_t massOmegaBckCandidate = 0.;
2141 Double_t widthOmegaBckCandidate = 0.;
2142
2143 omegaBckCandidate->GetMass(massOmegaBckCandidate,widthOmegaBckCandidate);
ef2e2748 2144
2145
48682642 2146 fHistograms->FillHistogram("ESD_Omega_Bck_InvMass_vs_Pt",massOmegaBckCandidate ,omegaBckCandidate->GetPt());
2147 fHistograms->FillHistogram("ESD_Omega_Bck_InvMass",massOmegaBckCandidate);
2148
2149 delete omegaBckCandidate;
2150 }
2151 }
2152 }
2153 } // end of checking if background calculation is available
6272370b 2154}
d7d7e825 2155
04bf4381 2156
77ac6f3e 2157void AliAnalysisTaskGammaConversion::AddOmegaToAOD(const AliKFParticle * const omegakf, Double_t mass, Int_t omegaDaughter, Int_t gammaDaughter) {
04bf4381 2158 //See header file for documentation
2159 AliGammaConversionAODObject omega;
77ac6f3e 2160 omega.SetPx(omegakf->GetPx());
2161 omega.SetPy(omegakf->GetPy());
2162 omega.SetPz(omegakf->GetPz());
04bf4381 2163 omega.SetChi2(omegakf->GetChi2());
77ac6f3e 2164 omega.SetE(omegakf->GetE());
04bf4381 2165 omega.SetIMass(mass);
2166 omega.SetLabel1(omegaDaughter);
6441e967 2167 // //dynamic_cast<AliAODPWG4Particle*>(fAODBranch->At(daughter1))->SetTagged(kTRUE);
04bf4381 2168 omega.SetLabel2(gammaDaughter);
6441e967 2169 AddToAODBranch(fAODOmega, omega);
04bf4381 2170}
2171
2172
2173
d7d7e825 2174void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
2175 // see header file for documentation
2176
6c84d371 2177 // for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){
2178 // for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){
037dc2db 2179
5ce758b0 2180 fESDEvent = fV0Reader->GetESDEvent();
2181
037dc2db 2182 if(fKFReconstructedGammasTClone->GetEntriesFast()>fV0Reader->GetNumberOfV0s()){
2183 cout<<"Warning, number of entries in the tclone is bigger than number of v0s"<<endl;
2184 }
2185
6c84d371 2186 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2187 for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();secondGammaIndex++){
d7d7e825 2188
6c84d371 2189 // AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];
2190 // AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];
a0b94e5c 2191
6c84d371 2192 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2193 AliKFParticle * twoGammaDecayCandidateDaughter1 = (AliKFParticle *)fKFReconstructedGammasTClone->At(secondGammaIndex);
a0b94e5c 2194
d7d7e825 2195 if(fElectronv1[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv1[firstGammaIndex]==fElectronv2[secondGammaIndex]){
2196 continue;
2197 }
2198 if(fElectronv2[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv2[firstGammaIndex]==fElectronv2[secondGammaIndex]){
2199 continue;
2200 }
a0b94e5c 2201
d7d7e825 2202 AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
2203
2204 Double_t massTwoGammaCandidate = 0.;
2205 Double_t widthTwoGammaCandidate = 0.;
2206 Double_t chi2TwoGammaCandidate =10000.;
2207 twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);
d707e3cf 2208 // if(twoGammaCandidate->GetNDF()>0){
2209 // chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();
2210 chi2TwoGammaCandidate = twoGammaCandidate->GetChi2();
d7d7e825 2211
70ef88b5 2212 fHistograms->FillHistogram("ESD_Mother_Chi2",chi2TwoGammaCandidate);
2e2da371 2213 if((chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
d7d7e825 2214
2215 TVector3 momentumVectorTwoGammaCandidate(twoGammaCandidate->GetPx(),twoGammaCandidate->GetPy(),twoGammaCandidate->GetPz());
2216 TVector3 spaceVectorTwoGammaCandidate(twoGammaCandidate->GetX(),twoGammaCandidate->GetY(),twoGammaCandidate->GetZ());
2217
2218 Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);
2219 Double_t rapidity;
96ade8ef 2220 if(twoGammaCandidate->GetE() - twoGammaCandidate->GetPz() <= 0 || twoGammaCandidate->GetE() + twoGammaCandidate->GetPz() <= 0){
2221 cout << "Error: |Pz| > E !!!! " << endl;
d7d7e825 2222 rapidity=0;
2223 }
2224 else{
2225 rapidity = 0.5*(TMath::Log((twoGammaCandidate->GetE() +twoGammaCandidate->GetPz()) / (twoGammaCandidate->GetE()-twoGammaCandidate->GetPz())));
2226 }
dc2883e4 2227
2228 if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut()){
2229 delete twoGammaCandidate;
2230 continue; // rapidity cut
2231 }
2232
2233
48682642 2234 Double_t alfa=0.0;
2235 if( (twoGammaDecayCandidateDaughter0->GetE()+twoGammaDecayCandidateDaughter1->GetE()) != 0){
2236 alfa=TMath::Abs((twoGammaDecayCandidateDaughter0->GetE()-twoGammaDecayCandidateDaughter1->GetE())
2237 /(twoGammaDecayCandidateDaughter0->GetE()+twoGammaDecayCandidateDaughter1->GetE()));
2238 }
d7d7e825 2239
1e7846f4 2240 if(openingAngleTwoGammaCandidate < fMinOpeningAngleGhostCut){
2241 delete twoGammaCandidate;
2242 continue; // minimum opening angle to avoid using ghosttracks
2243 }
2244
67381a40 2245 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
6de3471d 2246 fHistograms->FillHistogram("ESD_Mother_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);
2247 fHistograms->FillHistogram("ESD_Mother_Energy", twoGammaCandidate->GetE());
2248 fHistograms->FillHistogram("ESD_Mother_Pt", momentumVectorTwoGammaCandidate.Pt());
2249 fHistograms->FillHistogram("ESD_Mother_Eta", momentumVectorTwoGammaCandidate.Eta());
2250 fHistograms->FillHistogram("ESD_Mother_Rapid", rapidity);
2251 fHistograms->FillHistogram("ESD_Mother_Phi", spaceVectorTwoGammaCandidate.Phi());
2252 fHistograms->FillHistogram("ESD_Mother_Mass", massTwoGammaCandidate);
2253 fHistograms->FillHistogram("ESD_Mother_alfa", alfa);
9da64c46 2254 if(massTwoGammaCandidate>0.1 && massTwoGammaCandidate<0.15){
2255 fHistograms->FillHistogram("ESD_Mother_alfa_Pi0", alfa);
2256 }
6de3471d 2257 fHistograms->FillHistogram("ESD_Mother_R", spaceVectorTwoGammaCandidate.Pt()); // Pt in Space == R!!!
2258 fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), spaceVectorTwoGammaCandidate.Pt());
2259 fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());
2260 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2261 fHistograms->FillHistogram("ESD_Mother_InvMass",massTwoGammaCandidate);
9c1cb6f7 2262 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2263 }
3c45d101 2264 if(alfa<0.1){
2265 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_E_alpha",massTwoGammaCandidate ,twoGammaCandidate->GetE());
2266 }
6de3471d 2267
9640a3d1 2268
9c1cb6f7 2269 if(fCalculateBackground){
5ce758b0 2270 /* Kenneth, just for testing*/
2271 AliGammaConversionBGHandler * bgHandlerTest = fV0Reader->GetBGHandler();
2272
2273 Int_t zbin= bgHandlerTest->GetZBinIndex(fV0Reader->GetVertexZ());
2274 Int_t mbin=0;
2275 Int_t multKAA=0;
2276 if(fUseTrackMultiplicityForBG == kTRUE){
2277 multKAA=fV0Reader->CountESDTracks();
2278 mbin = bgHandlerTest->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2279 }
2280 else{// means we use #v0s for multiplicity
2281 multKAA=fV0Reader->GetNGoodV0s();
2282 mbin = bgHandlerTest->GetMultiplicityBinIndex(fV0Reader->GetNGoodV0s());
2283 }
2284 // cout<<"Filling bin number "<<zbin<<" and "<<mbin<<endl;
2285 // cout<<"Corresponding to z = "<<fV0Reader->GetVertexZ()<<" and m = "<<multKAA<<endl;
6de3471d 2286 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2287 fHistograms->FillHistogram(Form("%d%dESD_Mother_InvMass",zbin,mbin),massTwoGammaCandidate);
2288 fHistograms->FillHistogram(Form("%d%dESD_Mother_InvMass_vs_Pt",zbin,mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2289 /* end Kenneth, just for testing*/
2290 fHistograms->FillHistogram(Form("%dESD_Mother_InvMass_vs_Pt",mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2291 }
5ce758b0 2292 }
2293 /* if(fCalculateBackground){
9c1cb6f7 2294 AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
2295 Int_t mbin= bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2296 fHistograms->FillHistogram(Form("%dESD_Mother_InvMass_vs_Pt",mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
5ce758b0 2297 }*/
ebcfaa7e 2298 // if(fDoNeutralMesonV0MCCheck){
2299 if(fDoMCTruth){
037dc2db 2300 //Kenneth: Checking the eta of the gamma to check the difference between 0.9 and 1.2
2301 Int_t indexKF1 = fKFReconstructedGammasV0Index.at(firstGammaIndex);
2302 if(indexKF1<fV0Reader->GetNumberOfV0s()){
2303 fV0Reader->GetV0(indexKF1);//updates to the correct v0
2304 Double_t eta1 = fV0Reader->GetMotherCandidateEta();
2305 Bool_t isRealPi0=kFALSE;
10e3319b 2306 Bool_t isRealEta=kFALSE;
037dc2db 2307 Int_t gamma1MotherLabel=-1;
2308 if(fV0Reader->HasSameMCMother() == kTRUE){
2309 //cout<<"This v0 is a real v0!!!!"<<endl;
2310 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
2311 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
2312 if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
2313 if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){
2314 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
2315 gamma1MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
2316 }
2317 }
2318 }
2319 }
2320 Int_t indexKF2 = fKFReconstructedGammasV0Index.at(secondGammaIndex);
2321 if(indexKF1 == indexKF2){
2322 cout<<"index of the two KF particles are the same.... should not happen"<<endl;
2323 }
2324 if(indexKF2<fV0Reader->GetNumberOfV0s()){
2325 fV0Reader->GetV0(indexKF2);
2326 Double_t eta2 = fV0Reader->GetMotherCandidateEta();
2327 Int_t gamma2MotherLabel=-1;
2328 if(fV0Reader->HasSameMCMother() == kTRUE){
2329 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
2330 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
2331 if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
2332 if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){
2333 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
2334 gamma2MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
2335 }
2336 }
2337 }
2338 }
2339 if(gamma1MotherLabel>=0 && gamma1MotherLabel==gamma2MotherLabel){
2340 if(fV0Reader->CheckIfPi0IsMother(gamma1MotherLabel)){
2341 isRealPi0=kTRUE;
2342 }
10e3319b 2343 if(fV0Reader->CheckIfEtaIsMother(gamma1MotherLabel)){
2344 isRealEta=kTRUE;
2345 }
2346
037dc2db 2347 }
6de3471d 2348 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2349 if(TMath::Abs(eta1)>0.9 && TMath::Abs(eta2)>0.9){
2350 // fHistograms->FillHistogram("ESD_Mother_InvMass_1212",massTwoGammaCandidate);
2351 // fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2352 if(isRealPi0 || isRealEta){
2353 fHistograms->FillHistogram("ESD_TruePi0_InvMass_1212",massTwoGammaCandidate);
2354 fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_1212",openingAngleTwoGammaCandidate);
2355 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2356 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2357 fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
d707e3cf 2358 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2359 }
6de3471d 2360
2361 if(!isRealPi0 && !isRealEta){
2362 if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
2363 fHistograms->FillHistogram("ESD_TrueBckGG_InvMass",massTwoGammaCandidate);
2364 }else{
2365 fHistograms->FillHistogram("ESD_TrueBckCont_InvMass",massTwoGammaCandidate);
2366 }
bd6d9fa3 2367 }
bd6d9fa3 2368
6de3471d 2369 }
2370 else if(TMath::Abs(eta1)>0.9 || TMath::Abs(eta2)>0.9){
ebcfaa7e 2371 // fHistograms->FillHistogram("ESD_Mother_InvMass_0912",massTwoGammaCandidate);
2372 // fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
6de3471d 2373
2374 if(isRealPi0 || isRealEta){
2375 fHistograms->FillHistogram("ESD_TruePi0_InvMass_0912",massTwoGammaCandidate);
2376 fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_0912",openingAngleTwoGammaCandidate);
2377 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2378 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2379 fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
d707e3cf 2380 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2381 }
6de3471d 2382 if(!isRealPi0 && !isRealEta){
2383 if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
2384 fHistograms->FillHistogram("ESD_TrueBckGG_InvMass",massTwoGammaCandidate);
2385 }else{
2386 fHistograms->FillHistogram("ESD_TrueBckCont_InvMass",massTwoGammaCandidate);
2387 }
bd6d9fa3 2388 }
2389 }
6de3471d 2390 else{
ebcfaa7e 2391 // fHistograms->FillHistogram("ESD_Mother_InvMass_0909",massTwoGammaCandidate);
2392 // fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
6de3471d 2393 if(isRealPi0 || isRealEta){
2394 fHistograms->FillHistogram("ESD_TruePi0_InvMass_0909",massTwoGammaCandidate);
2395 fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_0909",openingAngleTwoGammaCandidate);
2396 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2397 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2398 fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
d707e3cf 2399 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
6de3471d 2400 if(gamma1MotherLabel > fV0Reader->GetMCStack()->GetNprimary()){
2401 fHistograms->FillHistogram("ESD_TruePi0Sec_InvMass",massTwoGammaCandidate);
2402 }
d707e3cf 2403 }
6de3471d 2404 if(!isRealPi0 && !isRealEta){
2405 if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
2406 fHistograms->FillHistogram("ESD_TrueBckGG_InvMass",massTwoGammaCandidate);
2407 }else{
2408 fHistograms->FillHistogram("ESD_TrueBckCont_InvMass",massTwoGammaCandidate);
2409 }
bd6d9fa3 2410 }
037dc2db 2411 }
2412 }
2413 }
2414 }
2415 }
6de3471d 2416 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
2417 if ( TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())<0.9 && TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())<0.9 ){
2418 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_Fiducial",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
2419 fHistograms->FillHistogram("ESD_Mother_InvMass_Fiducial",massTwoGammaCandidate);
2420 }
2421
2422 if(TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())>0.9 && TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())>0.9){
2423 fHistograms->FillHistogram("ESD_Mother_InvMass_1212",massTwoGammaCandidate);
2424 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2425 }
2426 else if(TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())>0.9 || TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())>0.9){
2427 fHistograms->FillHistogram("ESD_Mother_InvMass_0912",massTwoGammaCandidate);
2428 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2429 }
2430 else{
2431 fHistograms->FillHistogram("ESD_Mother_InvMass_0909",massTwoGammaCandidate);
2432 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
2433 }
ebcfaa7e 2434
6de3471d 2435 Double_t lowMassPi0=0.1;
2436 Double_t highMassPi0=0.15;
63477d82 2437 if ( fKFCreateAOD && (massTwoGammaCandidate > lowMassPi0) && (massTwoGammaCandidate < highMassPi0) ){
6de3471d 2438 new((*fKFReconstructedPi0sTClone)[fKFReconstructedPi0sTClone->GetEntriesFast()]) AliKFParticle(*twoGammaCandidate);
2439 fGammav1.push_back(firstGammaIndex);
2440 fGammav2.push_back(secondGammaIndex);
22eae5f9 2441 AddPionToAOD(twoGammaCandidate, massTwoGammaCandidate, firstGammaIndex, secondGammaIndex);
6de3471d 2442 }
ebcfaa7e 2443 }
6272370b 2444
2e2da371 2445 }
d707e3cf 2446 //}
6de3471d 2447 delete twoGammaCandidate;
d7d7e825 2448 }
2449 }
2450}
2451
22eae5f9 2452void AliAnalysisTaskGammaConversion::AddPionToAOD(AliKFParticle * pionkf, Double_t mass, Int_t daughter1, Int_t daughter2) {
04bf4381 2453 //See header file for documentation
22eae5f9 2454 if(fOutputAODClassName.Contains("AODObject")) {
2455 AliGammaConversionAODObject pion;
2456 pion.SetPx(pionkf->GetPx());
2457 pion.SetPy(pionkf->GetPy());
2458 pion.SetPz(pionkf->GetPz());
2459 pion.SetChi2(pionkf->GetChi2());
2460 pion.SetE(pionkf->GetE());
2461 pion.SetIMass(mass);
2462 pion.SetLabel1(daughter1);
2463 pion.SetLabel2(daughter2);
2464 AddToAODBranch(fAODPi0, pion);
2465 } else {
2466 TLorentzVector momentum(pionkf->Px(),pionkf->Py(),pionkf->Pz(), pionkf->E());
2467 AliAODConversionParticle pion = AliAODConversionParticle(momentum);
2468 pion.SetTrackLabels( daughter1, daughter2 );
2469 pion.SetChi2(pionkf->GetChi2());
2470 AddToAODBranch(fAODPi0, pion);
2471
2472 }
04bf4381 2473}
2474
2475
48682642 2476/*
77ac6f3e 2477 void AliAnalysisTaskGammaConversion::ProcessConvPHOSGammasForNeutralMesonAnalysis(){
2478
6272370b 2479 // see header file for documentation
2480 // Analyse Pi0 with one photon from Phos and 1 photon from conversions
2481
2482
2483
2484 Double_t vtx[3];
2485 vtx[0] = fV0Reader->GetPrimaryVertex()->GetX();
2486 vtx[1] = fV0Reader->GetPrimaryVertex()->GetY();
2487 vtx[2] = fV0Reader->GetPrimaryVertex()->GetZ();
2488
2489
2490 // Loop over all CaloClusters and consider only the PHOS ones:
2491 AliESDCaloCluster *clu;
2492 TLorentzVector pPHOS;
2493 TLorentzVector gammaPHOS;
2494 TLorentzVector gammaGammaConv;
2495 TLorentzVector pi0GammaConvPHOS;
2496 TLorentzVector gammaGammaConvBck;
2497 TLorentzVector pi0GammaConvPHOSBck;
2498
2499
2500 for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
77ac6f3e 2501 clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
2502 if ( !clu->IsPHOS() || clu->E()<0.1 ) continue;
2503 clu ->GetMomentum(pPHOS ,vtx);
2504 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2505 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2506 gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz(),0.);
2507 gammaPHOS.SetXYZM(pPHOS.Px(),pPHOS.Py(),pPHOS.Pz(),0.);
2508 pi0GammaConvPHOS=gammaGammaConv+gammaPHOS;
2509 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS",pi0GammaConvPHOS.M());
2510 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvPHOS",pi0GammaConvPHOS.M(),pi0GammaConvPHOS.Pt());
6272370b 2511
77ac6f3e 2512 TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz());
2513 TVector3 v3D1(gammaPHOS.Px(),gammaPHOS.Py(),gammaPHOS.Pz());
2514 Double_t opanConvPHOS= v3D0.Angle(v3D1);
2515 if ( opanConvPHOS < 0.35){
2516 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanLow",pi0GammaConvPHOS.M());
2517 }else{
2518 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanHigh",pi0GammaConvPHOS.M());
2519 }
6272370b 2520
77ac6f3e 2521 }
2522
2523 // Now the LorentVector pPHOS is obtained and can be paired with the converted proton
6272370b 2524 }
2525 //==== End of the PHOS cluster selection ============
2526 TLorentzVector pEMCAL;
2527 TLorentzVector gammaEMCAL;
2528 TLorentzVector pi0GammaConvEMCAL;
2529 TLorentzVector pi0GammaConvEMCALBck;
2530
77ac6f3e 2531 for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
2532 clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
2533 if ( !clu->IsEMCAL() || clu->E()<0.1 ) continue;
2534 if (clu->GetNCells() <= 1) continue;
2535 if ( clu->GetTOF()*1e9 < 550 || clu->GetTOF()*1e9 > 750) continue;
6272370b 2536
77ac6f3e 2537 clu ->GetMomentum(pEMCAL ,vtx);
2538 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2539 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
2540 gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),
2541 twoGammaDecayCandidateDaughter0->Py(),
2542 twoGammaDecayCandidateDaughter0->Pz(),0.);
2543 gammaEMCAL.SetXYZM(pEMCAL.Px(),pEMCAL.Py(),pEMCAL.Pz(),0.);
2544 pi0GammaConvEMCAL=gammaGammaConv+gammaEMCAL;
2545 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL",pi0GammaConvEMCAL.M());
2546 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL",pi0GammaConvEMCAL.M(),pi0GammaConvEMCAL.Pt());
2547 TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),
2548 twoGammaDecayCandidateDaughter0->Py(),
2549 twoGammaDecayCandidateDaughter0->Pz());
2550 TVector3 v3D1(gammaEMCAL.Px(),gammaEMCAL.Py(),gammaEMCAL.Pz());
6272370b 2551
77ac6f3e 2552
2553 Double_t opanConvEMCAL= v3D0.Angle(v3D1);
2554 if ( opanConvEMCAL < 0.35){
2555 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanLow",pi0GammaConvEMCAL.M());
2556 }else{
2557 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanHigh",pi0GammaConvEMCAL.M());
2558 }
2559
2560 }
2561 if(fCalculateBackground){
2562 for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2563 AliGammaConversionKFVector * previousEventV0s = fV0Reader->GetBGGoodV0s(nEventsInBG);
2564 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2565 AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2566 gammaGammaConvBck.SetXYZM(previousGoodV0.Px(),
2567 previousGoodV0.Py(),
2568 previousGoodV0.Pz(),0.);
2569 pi0GammaConvEMCALBck=gammaGammaConvBck+gammaEMCAL;
2570 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M());
2571 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M(),
2572 pi0GammaConvEMCALBck.Pt());
2573 }
2574 }
48682642 2575
77ac6f3e 2576 // Now the LorentVector pEMCAL is obtained and can be paired with the converted proton
2577 } // end of checking if background photons are available
2578 }
6272370b 2579 //==== End of the PHOS cluster selection ============
77ac6f3e 2580
2581 }
48682642 2582*/
6272370b 2583
77ac6f3e 2584void AliAnalysisTaskGammaConversion::MoveParticleAccordingToVertex(AliKFParticle * particle,const AliGammaConversionBGHandler::GammaConversionVertex *vertex){
5ce758b0 2585 //see header file for documentation
2586
77ac6f3e 2587 Double_t dx = vertex->fX - fESDEvent->GetPrimaryVertex()->GetX();
2588 Double_t dy = vertex->fY - fESDEvent->GetPrimaryVertex()->GetY();
2589 Double_t dz = vertex->fZ - fESDEvent->GetPrimaryVertex()->GetZ();
5ce758b0 2590
2591 // cout<<"dx, dy, dz: ["<<dx<<","<<dy<<","<<dz<<"]"<<endl;
77ac6f3e 2592 particle->X() = particle->GetX() - dx;
2593 particle->Y() = particle->GetY() - dy;
2594 particle->Z() = particle->GetZ() - dz;
5ce758b0 2595}
2596
111d75df 2597void AliAnalysisTaskGammaConversion::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle){
77ac6f3e 2598 // Rotate the kf particle
111d75df 2599 Double_t c = cos(angle);
2600 Double_t s = sin(angle);
2601
77ac6f3e 2602 Double_t mA[7][ 7];
111d75df 2603 for( Int_t i=0; i<7; i++ ){
2604 for( Int_t j=0; j<7; j++){
77ac6f3e 2605 mA[i][j] = 0;
111d75df 2606 }
2607 }
2608 for( int i=0; i<7; i++ ){
77ac6f3e 2609 mA[i][i] = 1;
111d75df 2610 }
77ac6f3e 2611 mA[0][0] = c; mA[0][1] = s;
2612 mA[1][0] = -s; mA[1][1] = c;
2613 mA[3][3] = c; mA[3][4] = s;
2614 mA[4][3] = -s; mA[4][4] = c;
111d75df 2615
77ac6f3e 2616 Double_t mAC[7][7];
2617 Double_t mAp[7];
111d75df 2618
2619 for( Int_t i=0; i<7; i++ ){
77ac6f3e 2620 mAp[i] = 0;
111d75df 2621 for( Int_t k=0; k<7; k++){
77ac6f3e 2622 mAp[i]+=mA[i][k] * kfParticle->GetParameter(k);
111d75df 2623 }
2624 }
2625
2626 for( Int_t i=0; i<7; i++){
77ac6f3e 2627 kfParticle->Parameter(i) = mAp[i];
111d75df 2628 }
2629
2630 for( Int_t i=0; i<7; i++ ){
2631 for( Int_t j=0; j<7; j++ ){
77ac6f3e 2632 mAC[i][j] = 0;
111d75df 2633 for( Int_t k=0; k<7; k++ ){
77ac6f3e 2634 mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j);
111d75df 2635 }
2636 }
2637 }
2638
2639 for( Int_t i=0; i<7; i++ ){
2640 for( Int_t j=0; j<=i; j++ ){
2641 Double_t xx = 0;
2642 for( Int_t k=0; k<7; k++){
77ac6f3e 2643 xx+= mAC[i][k]*mA[j][k];
111d75df 2644 }
2645 kfParticle->Covariance(i,j) = xx;
2646 }
2647 }
2648}
2649
2650
d7d7e825 2651void AliAnalysisTaskGammaConversion::CalculateBackground(){
2652 // see header file for documentation
5e55d806 2653
2654
2655 TClonesArray * currentEventV0s = fV0Reader->GetCurrentEventGoodV0s();
2656
5ce758b0 2657 AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
2658
2659 Int_t zbin= bgHandler->GetZBinIndex(fV0Reader->GetVertexZ());
2660 Int_t mbin = 0;
2661 if(fUseTrackMultiplicityForBG == kTRUE){
2662 mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2663 }
2664 else{
2665 mbin = bgHandler->GetMultiplicityBinIndex(fV0Reader->GetNGoodV0s());
2666 }
2667
111d75df 2668 if(fDoRotation == kTRUE){
3c45d101 2669 TRandom3 *random = new TRandom3();
111d75df 2670 for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2671 AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
2672 for(Int_t iCurrent2=iCurrent+1;iCurrent2<currentEventV0s->GetEntriesFast();iCurrent2++){
77ac6f3e 2673 for(Int_t nRandom=0;nRandom<fNRandomEventsForBG;nRandom++){
111d75df 2674
2675 AliKFParticle currentEventGoodV02 = *(AliKFParticle *)(currentEventV0s->At(iCurrent2));
3bfbe89a 2676
2677 if(fCheckBGProbability == kTRUE){
2678 Double_t massBGprob =0.;
2679 Double_t widthBGprob = 0.;
2680 AliKFParticle *backgroundCandidateProb = new AliKFParticle(currentEventGoodV0,currentEventGoodV02);
2681 backgroundCandidateProb->GetMass(massBGprob,widthBGprob);
2682 if(massBGprob>0.1 && massBGprob<0.14){
2683 if(random->Rndm()>bgHandler->GetBGProb(zbin,mbin)){
2684 delete backgroundCandidateProb;
2685 continue;
2686 }
2687 }
2688 delete backgroundCandidateProb;
2689 }
111d75df 2690
77ac6f3e 2691 Double_t nRadiansPM = fNDegreesPMBackground*TMath::Pi()/180;
111d75df 2692
2693 Double_t rotationValue = random->Rndm()*2*nRadiansPM + TMath::Pi()-nRadiansPM;
5ce758b0 2694
111d75df 2695 RotateKFParticle(&currentEventGoodV02,rotationValue);
5ce758b0 2696
111d75df 2697 AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,currentEventGoodV02);
5e55d806 2698
5ce758b0 2699 Double_t massBG =0.;
2700 Double_t widthBG = 0.;
2701 Double_t chi2BG =10000.;
2702 backgroundCandidate->GetMass(massBG,widthBG);
111d75df 2703
111d75df 2704 // if(backgroundCandidate->GetNDF()>0){
2705 chi2BG = backgroundCandidate->GetChi2();
2706 if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
2707
2708 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2709 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2710
2711 Double_t openingAngleBG = currentEventGoodV0.GetAngle(currentEventGoodV02);
2712
2713 Double_t rapidity;
2714 if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
2715 else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2716
dc2883e4 2717 if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ){
2718 delete backgroundCandidate;
2719 continue; // rapidity cut
2720 }
2721
111d75df 2722
2723 Double_t alfa=0.0;
2724 if( (currentEventGoodV0.GetE()+currentEventGoodV02.GetE()) != 0){
2725 alfa=TMath::Abs((currentEventGoodV0.GetE()-currentEventGoodV02.GetE())
2726 /(currentEventGoodV0.GetE()+currentEventGoodV02.GetE()));
2727 }
2728
2729
2730 if(openingAngleBG < fMinOpeningAngleGhostCut ){
2731 delete backgroundCandidate;
2732 continue; // minimum opening angle to avoid using ghosttracks
2733 }
2734
2735 // original
111d75df 2736 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
6de3471d 2737 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2738 fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2739 fHistograms->FillHistogram("ESD_Background_Pt", momentumVectorbackgroundCandidate.Pt());
2740 fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2741 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2742 fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2743 fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2744 fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2745 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2746 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2747 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2748 fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
111d75df 2749 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
6de3471d 2750
2751
2752 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(currentEventGoodV02.GetEta())<0.9 ){
2753 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
2754 fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
2755 }
2756
2757 fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
2758 fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
2759 fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin), momentumVectorbackgroundCandidate.Pt());
2760 fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
2761 fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
2762 fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
2763 fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
2764 fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2765 fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2766 fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
2767 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2768 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
2769
2770 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(currentEventGoodV02.GetEta())<0.9 ){
2771 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2772 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
2773 }
111d75df 2774 }
3c45d101 2775 if(alfa<0.1){
2776 fHistograms->FillHistogram("ESD_Background_InvMass_vs_E_alpha",massBG ,backgroundCandidate->GetE());
2777 }
2778
2e2da371 2779 }
111d75df 2780 //}
5ce758b0 2781 delete backgroundCandidate;
2782 }
2783 }
2784 }
111d75df 2785 delete random;
5ce758b0 2786 }
111d75df 2787 else{ // means no rotation
2788 AliGammaConversionBGHandler::GammaConversionVertex *bgEventVertex = NULL;
2789
2790 if(fUseTrackMultiplicityForBG){
2791 // cout<<"Using charged track multiplicity for background calculation"<<endl;
2792 for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
5ce758b0 2793
111d75df 2794 AliGammaConversionKFVector * previousEventV0s = bgHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);//fV0Reader->GetBGGoodV0s(nEventsInBG);
2795
5ce758b0 2796 if(fMoveParticleAccordingToVertex == kTRUE){
2797 bgEventVertex = bgHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
2798 }
2799
2800 for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2801 AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
2802 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2803 AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
111d75df 2804 AliKFParticle previousGoodV0test = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
5ce758b0 2805
111d75df 2806 //cout<<"Primary Vertex event: ["<<fESDEvent->GetPrimaryVertex()->GetX()<<","<<fESDEvent->GetPrimaryVertex()->GetY()<<","<<fESDEvent->GetPrimaryVertex()->GetZ()<<"]"<<endl;
2807 //cout<<"BG prim Vertex event: ["<<bgEventVertex->fX<<","<<bgEventVertex->fY<<","<<bgEventVertex->fZ<<"]"<<endl;
2808
2809 //cout<<"XYZ of particle before transport: ["<<previousGoodV0.X()<<","<<previousGoodV0.Y()<<","<<previousGoodV0.Z()<<"]"<<endl;
5ce758b0 2810 if(fMoveParticleAccordingToVertex == kTRUE){
2811 MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
2812 }
111d75df 2813 //cout<<"XYZ of particle after transport: ["<<previousGoodV0.X()<<","<<previousGoodV0.Y()<<","<<previousGoodV0.Z()<<"]"<<endl;
5ce758b0 2814
2815 AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
111d75df 2816
5ce758b0 2817 Double_t massBG =0.;
2818 Double_t widthBG = 0.;
2819 Double_t chi2BG =10000.;
2820 backgroundCandidate->GetMass(massBG,widthBG);
111d75df 2821 // if(backgroundCandidate->GetNDF()>0){
2822 // chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
2823 chi2BG = backgroundCandidate->GetChi2();
2824 if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
2825
2826 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2827 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2828
2829 Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2830
2831 Double_t rapidity;
2832
2833 if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() <= 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() <= 0){
2834 cout << "Error: |Pz| > E !!!! " << endl;
2835 rapidity=0;
2836 } else {
2837 rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2838 }
dc2883e4 2839 if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ){
2840 delete backgroundCandidate;
2841 continue; // rapidity cut
2842 }
2843
2844
111d75df 2845 Double_t alfa=0.0;
2846 if( (currentEventGoodV0.GetE()+previousGoodV0.GetE()) != 0){
2847 alfa=TMath::Abs((currentEventGoodV0.GetE()-previousGoodV0.GetE())
2848 /(currentEventGoodV0.GetE()+previousGoodV0.GetE()));
2849 }
2850
2851
2852 if(openingAngleBG < fMinOpeningAngleGhostCut ){
2853 delete backgroundCandidate;
2854 continue; // minimum opening angle to avoid using ghosttracks
2855 }
2856
2857 // original
111d75df 2858 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
6de3471d 2859 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2860 fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2861 fHistograms->FillHistogram("ESD_Background_Pt", momentumVectorbackgroundCandidate.Pt());
2862 fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2863 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2864 fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2865 fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2866 fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2867 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2868 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2869 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2870 fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
111d75df 2871 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
6de3471d 2872
2873
2874 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2875 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
2876 fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
2877 }
2878
2879 // test
2880 fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
2881 fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
2882 fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin), momentumVectorbackgroundCandidate.Pt());
2883 fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
2884 fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
2885 fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
2886 fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
2887 fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2888 fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2889 fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
2890 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2891 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
2892
2893 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2894 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2895 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
2896 }
2897 // }
111d75df 2898 }
3c45d101 2899 if(alfa<0.1){
2900 fHistograms->FillHistogram("ESD_Background_InvMass_vs_E_alpha",massBG ,backgroundCandidate->GetE());
2901 }
2902
111d75df 2903 }
2904 delete backgroundCandidate;
2905 }
2906 }
2907 }
2908 }
2909 else{ // means using #V0s for multiplicity
2910
2911 // cout<<"Using the v0 multiplicity to calculate background"<<endl;
2912
2913 fHistograms->FillHistogram("ESD_Background_z_m",zbin,mbin);
2914 fHistograms->FillHistogram("ESD_Mother_multpilicityVSv0s",fV0Reader->CountESDTracks(),fV0Reader->GetNumberOfV0s());
2915
2916 for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2917 AliGammaConversionKFVector * previousEventV0s = bgHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);// fV0Reader->GetBGGoodV0s(nEventsInBG);
2918 if(previousEventV0s){
2919
2920 if(fMoveParticleAccordingToVertex == kTRUE){
2921 bgEventVertex = bgHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
2922 }
2923
2924 for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2925 AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
2926 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2927 AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2928
2929 if(fMoveParticleAccordingToVertex == kTRUE){
2930 MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
5ce758b0 2931 }
111d75df 2932
2933 AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
2934 Double_t massBG =0.;
2935 Double_t widthBG = 0.;
2936 Double_t chi2BG =10000.;
2937 backgroundCandidate->GetMass(massBG,widthBG);
2938 /* if(backgroundCandidate->GetNDF()>0){
2939 chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
2940 {//remember to remove
2941 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2942 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2943
2944 Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2945 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle_nochi2", openingAngleBG);
2946 }
2947 */
2e2da371 2948 chi2BG = backgroundCandidate->GetChi2();
2949 if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()) || fApplyChi2Cut == kFALSE){
5ce758b0 2950 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2951 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2952
2953 Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2954
2955 Double_t rapidity;
2956 if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
2957 else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2958
dc2883e4 2959 if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ){
2960 delete backgroundCandidate;
2961 continue; // rapidity cut
2962 }
2963
2964
5ce758b0 2965 Double_t alfa=0.0;
2966 if( (currentEventGoodV0.GetE()+previousGoodV0.GetE()) != 0){
2967 alfa=TMath::Abs((currentEventGoodV0.GetE()-previousGoodV0.GetE())
2968 /(currentEventGoodV0.GetE()+previousGoodV0.GetE()));
2969 }
2970
2971
2972 if(openingAngleBG < fMinOpeningAngleGhostCut ){
2973 delete backgroundCandidate;
2974 continue; // minimum opening angle to avoid using ghosttracks
2975 }
2976
5ce758b0 2977 if(alfa>fV0Reader->GetAlphaMinCutMeson() && alfa<fV0Reader->GetAlphaCutMeson()){
6de3471d 2978 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2979 fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2980 fHistograms->FillHistogram("ESD_Background_Pt", momentumVectorbackgroundCandidate.Pt());
2981 fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2982 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2983 fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2984 fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2985 fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2986 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2987 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2988 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2989 fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
2990
2991
5ce758b0 2992 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
3c45d101 2993
6de3471d 2994
2995 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2996 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
2997 fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
2998 }
2999
3000 if(massBG>0.5 && massBG<0.6){
3001 fHistograms->FillHistogram("ESD_Background_alfa_pt0506",momentumVectorbackgroundCandidate.Pt(),alfa);
3002 }
3003 if(massBG>0.3 && massBG<0.4){
3004 fHistograms->FillHistogram("ESD_Background_alfa_pt0304",momentumVectorbackgroundCandidate.Pt(),alfa);
3005 }
3006
3007 // test
3008 fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
3009 fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
3010 fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin), momentumVectorbackgroundCandidate.Pt());
3011 fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
3012 fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
3013 fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
3014 fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
3015 fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
3016 fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
3017 fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
3018 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
3019 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
3020
3021 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
3022 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
3023 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
3024 }
5ce758b0 3025 }
5ce758b0 3026
6de3471d 3027 if(alfa<0.1){
3028 fHistograms->FillHistogram("ESD_Background_InvMass_vs_E_alpha",massBG ,backgroundCandidate->GetE());
5ce758b0 3029 }
2e2da371 3030 // }
5ce758b0 3031 }
111d75df 3032 delete backgroundCandidate;
3033 }
5ce758b0 3034 }
3035 }
d7d7e825 3036 }
111d75df 3037 } // end else (means use #v0s as multiplicity)
3038 } // end no rotation
d7d7e825 3039}
3040
3041
d7d7e825 3042void AliAnalysisTaskGammaConversion::ProcessGammasForGammaJetAnalysis(){
3043 //ProcessGammasForGammaJetAnalysis
a0b94e5c 3044
d7d7e825 3045 Double_t distIsoMin;
a0b94e5c 3046
d7d7e825 3047 CreateListOfChargedParticles();
a0b94e5c 3048
3049
6c84d371 3050 // for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
3051 for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
3052 AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
01b7fdcc 3053 TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
a0b94e5c 3054
01b7fdcc 3055 if( momentumVectorCurrentGamma.Pt()> fMinPtForGammaJet){
d7d7e825 3056 distIsoMin=GetMinimumDistanceToCharge(gammaIndex);
a0b94e5c 3057
d7d7e825 3058 if (distIsoMin > fMinIsoConeSize && fLeadingChargedIndex>=0){
01b7fdcc 3059 CalculateJetCone(gammaIndex);
d7d7e825 3060 }
3061 }
3062 }
3063}
3064
6272370b 3065//____________________________________________________________________
77ac6f3e 3066Bool_t AliAnalysisTaskGammaConversion::IsGoodImpPar(const AliESDtrack *const track)
6272370b 3067{
3068//
3069// check whether particle has good DCAr(Pt) impact
3070// parameter. Only for TPC+ITS tracks (7*sigma cut)
3071// Origin: Andrea Dainese
3072//
3073
3074Float_t d0z0[2],covd0z0[3];
3075track->GetImpactParameters(d0z0,covd0z0);
3076Float_t sigma= 0.0050+0.0060/TMath::Power(track->Pt(),0.9);
3077Float_t d0max = 7.*sigma;
3078if(TMath::Abs(d0z0[0]) < d0max) return kTRUE;
3079
3080return kFALSE;
3081}
3082
3083
d7d7e825 3084void AliAnalysisTaskGammaConversion::CreateListOfChargedParticles(){
3085 // CreateListOfChargedParticles
a0b94e5c 3086
d7d7e825 3087 fESDEvent = fV0Reader->GetESDEvent();
037dc2db 3088 Int_t numberOfESDTracks=0;
d7d7e825 3089 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
3090 AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
a0b94e5c 3091
d7d7e825 3092 if(!curTrack){
3093 continue;
3094 }
d707e3cf 3095 // Not needed if Standard function used.
3096// if(!IsGoodImpPar(curTrack)){
3097// continue;
3098// }
a0b94e5c 3099
d7d7e825 3100 if(fEsdTrackCuts->AcceptTrack(curTrack) ){
6c84d371 3101 new((*fChargedParticles)[fChargedParticles->GetEntriesFast()]) AliESDtrack(*curTrack);
3102 // fChargedParticles.push_back(curTrack);
d7d7e825 3103 fChargedParticlesId.push_back(iTracks);
037dc2db 3104 numberOfESDTracks++;
d7d7e825 3105 }
3106 }
6746e1e1 3107// Moved to UserExec using CountAcceptedTracks function. runjet is not needed by default
3108// fHistograms->FillHistogram("ESD_NumberOfGoodESDTracks",numberOfESDTracks);
3109// cout<<"esdtracks::"<< numberOfESDTracks<<endl;
3110// if (fV0Reader->GetNumberOfContributorsVtx()>=1){
3111// fHistograms->FillHistogram("ESD_NumberOfGoodESDTracksVtx",numberOfESDTracks);
3112// }
d7d7e825 3113}
9c1cb6f7 3114void AliAnalysisTaskGammaConversion::RecalculateV0ForGamma(){
77ac6f3e 3115 //recalculates v0 for gamma
3116
9c1cb6f7 3117 Double_t massE=0.00051099892;
3118 TLorentzVector curElecPos;
3119 TLorentzVector curElecNeg;
3120 TLorentzVector curGamma;
3121
3122 TLorentzVector curGammaAt;
3123 TLorentzVector curElecPosAt;
3124 TLorentzVector curElecNegAt;
3125 AliKFVertex primVtxGamma(*(fESDEvent->GetPrimaryVertex()));
3126 AliKFVertex primVtxImprovedGamma = primVtxGamma;
3127
3128 const AliESDVertex *vtxT3D=fESDEvent->GetPrimaryVertex();
3129
3130 Double_t xPrimaryVertex=vtxT3D->GetXv();
3131 Double_t yPrimaryVertex=vtxT3D->GetYv();
3132 Double_t zPrimaryVertex=vtxT3D->GetZv();
70ef88b5 3133 // Float_t primvertex[3]={xPrimaryVertex,yPrimaryVertex,zPrimaryVertex};
9c1cb6f7 3134
3135 Float_t nsigmaTPCtrackPos;
3136 Float_t nsigmaTPCtrackNeg;
3137 Float_t nsigmaTPCtrackPosToPion;
3138 Float_t nsigmaTPCtrackNegToPion;
3139 AliKFParticle* negKF=NULL;
3140 AliKFParticle* posKF=NULL;
3141
3142 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
3143 AliESDtrack* posTrack = fESDEvent->GetTrack(iTracks);
3144 if(!posTrack){
3145 continue;
3146 }
3147 if (posKF) delete posKF; posKF=NULL;
3148 if(posTrack->GetSign()<0) continue;
3149 if(!(posTrack->GetStatus() & AliESDtrack::kTPCrefit))continue;
3150 if(posTrack->GetKinkIndex(0)>0 ) continue;
3151 if(posTrack->GetNcls(1)<50)continue;
3152 Double_t momPos[3];
3153 // posTrack->GetConstrainedPxPyPz(momPos);
3154 posTrack->GetPxPyPz(momPos);
3155 AliESDtrack *ptrk=fESDEvent->GetTrack(iTracks);
3156 curElecPos.SetXYZM(momPos[0],momPos[1],momPos[2],massE);
3157 if(TMath::Abs(curElecPos.Eta())<0.9) continue;
3158 posKF = new AliKFParticle( *(posTrack),-11);
3159
3160 nsigmaTPCtrackPos = fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kElectron);
3161 nsigmaTPCtrackPosToPion = fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kPion);
3162
3163 if ( nsigmaTPCtrackPos>5.|| nsigmaTPCtrackPos<-2.){
3164 continue;
3165 }
3166
3167 if(pow((momPos[0]*momPos[0]+momPos[1]*momPos[1]+momPos[2]*momPos[2]),0.5)>0.5 && nsigmaTPCtrackPosToPion<1){
3168 continue;
3169 }
3170
3171
3172
3173 for(Int_t jTracks = 0; jTracks < fESDEvent->GetNumberOfTracks(); jTracks++){
3174 AliESDtrack* negTrack = fESDEvent->GetTrack(jTracks);
3175 if(!negTrack){
3176 continue;
3177 }
3178 if (negKF) delete negKF; negKF=NULL;
3179 if(negTrack->GetSign()>0) continue;
3180 if(!(negTrack->GetStatus() & AliESDtrack::kTPCrefit))continue;
3181 if(negTrack->GetKinkIndex(0)>0 ) continue;
3182 if(negTrack->GetNcls(1)<50)continue;
3183 Double_t momNeg[3];
3184 // negTrack->GetConstrainedPxPyPz(momNeg);
3185 negTrack->GetPxPyPz(momNeg);
3186
3187 nsigmaTPCtrackNeg = fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kElectron);
3188 nsigmaTPCtrackNegToPion = fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kPion);
3189 if ( nsigmaTPCtrackNeg>5. || nsigmaTPCtrackNeg<-2.){
3190 continue;
3191 }
3192 if(pow((momNeg[0]*momNeg[0]+momNeg[1]*momNeg[1]+momNeg[2]*momNeg[2]),0.5)>0.5 && nsigmaTPCtrackNegToPion<1){
3193 continue;
3194 }
3195 AliESDtrack *ntrk=fESDEvent->GetTrack(jTracks);
3196 curElecNeg.SetXYZM(momNeg[0],momNeg[1],momNeg[2],massE);
3197 if(TMath::Abs(curElecNeg.Eta())<0.9) continue;
3198 negKF = new AliKFParticle( *(negTrack) ,11);
3199
3200 Double_t b=fESDEvent->GetMagneticField();
3201 Double_t xn, xp, dca=ntrk->GetDCA(ptrk,b,xn,xp);
3202 AliExternalTrackParam nt(*ntrk), pt(*ptrk);
3203 nt.PropagateTo(xn,b); pt.PropagateTo(xp,b);
3204
3205
3206 //--- Like in ITSV0Finder
3207 AliExternalTrackParam ntAt0(*ntrk), ptAt0(*ptrk);
3208 Double_t xxP,yyP,alphaP;
3209 Double_t rP[3];
3210
3211 // if (!ptAt0.GetGlobalXYZat(ptAt0->GetX(),xxP,yyP,zzP)) continue;
3212 if (!ptAt0.GetXYZAt(ptAt0.GetX(),b,rP)) continue;
3213 xxP=rP[0];
3214 yyP=rP[1];
3215 alphaP = TMath::ATan2(yyP,xxP);
3216
3217
3218 ptAt0.Propagate(alphaP,0,b);
3219 Float_t ptfacP = (1.+100.*TMath::Abs(ptAt0.GetC(b)));
3220
3221 // Double_t distP = ptAt0.GetY();
3222 Double_t normP = ptfacP*TMath::Sqrt(ptAt0.GetSigmaY2());
3223 Double_t normdist0P = TMath::Abs(ptAt0.GetY()/normP);
3224 Double_t normdist1P = TMath::Abs((ptAt0.GetZ()-zPrimaryVertex)/(ptfacP*TMath::Sqrt(ptAt0.GetSigmaZ2())));
3225 Double_t normdistP = TMath::Sqrt(normdist0P*normdist0P+normdist1P*normdist1P);
3226
3227
3228 Double_t xxN,yyN,alphaN;
3229 Double_t rN[3];
3230 // if (!ntAt0.GetGlobalXYZat(ntAt0->GetX(),xxN,yyN,zzN)) continue;
3231 if (!ntAt0.GetXYZAt(ntAt0.GetX(),b,rN)) continue;
3232 xxN=rN[0];
3233 yyN=rN[1];
3234
3235 alphaN = TMath::ATan2(yyN,xxN);
3236
3237 ntAt0.Propagate(alphaN,0,b);
3238
3239 Float_t ptfacN = (1.+100.*TMath::Abs(ntAt0.GetC(b)));
3240 // Double_t distN = ntAt0.GetY();
3241 Double_t normN = ptfacN*TMath::Sqrt(ntAt0.GetSigmaY2());
3242 Double_t normdist0N = TMath::Abs(ntAt0.GetY()/normN);
3243 Double_t normdist1N = TMath::Abs((ntAt0.GetZ()-zPrimaryVertex)/(ptfacN*TMath::Sqrt(ntAt0.GetSigmaZ2())));
3244 Double_t normdistN = TMath::Sqrt(normdist0N*normdist0N+normdist1N*normdist1N);
3245
3246 //-----------------------------
3247
3248 Double_t momNegAt[3];
3249 nt.GetPxPyPz(momNegAt);
3250 curElecNegAt.SetXYZM(momNegAt[0],momNegAt[1],momNegAt[2],massE);
3251
3252 Double_t momPosAt[3];
3253 pt.GetPxPyPz(momPosAt);
3254 curElecPosAt.SetXYZM(momPosAt[0],momPosAt[1],momPosAt[2],massE);
3255 if(dca>1){
3256 continue;
3257 }
3258
3259 // Double_t dneg= negTrack->GetD(xPrimaryVertex,yPrimaryVertex,b);
3260 // Double_t dpos= posTrack->GetD(xPrimaryVertex,yPrimaryVertex,b);
3261 AliESDv0 vertex(nt,jTracks,pt,iTracks);
3262
3263
3264 Float_t cpa=vertex.GetV0CosineOfPointingAngle(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex);
3265
9c1cb6f7 3266
70ef88b5 3267
9c1cb6f7 3268 // cout<< "v0Rr::"<< v0Rr<<endl;
3269 // if (pvertex.GetRr()<0.5){
3270 // continue;
3271 //}
3272 // cout<<"vertex.GetChi2V0()"<<vertex.GetChi2V0()<<endl;
3273 if(cpa<0.9)continue;
3274 // if (vertex.GetChi2V0() > 30) continue;
3275 // cout<<"xp+xn::"<<xp<<" "<<xn<<endl;
3276 if ((xn+xp) < 0.4) continue;
3277 if (TMath::Abs(ntrk->GetD(xPrimaryVertex,yPrimaryVertex,b))<0.05)
3278 if (TMath::Abs(ptrk->GetD(xPrimaryVertex,yPrimaryVertex,b))<0.05) continue;
3279
3280 //cout<<"pass"<<endl;
3281
3282 AliKFParticle v0GammaC;
3283 v0GammaC+=(*negKF);
3284 v0GammaC+=(*posKF);
3285 v0GammaC.SetMassConstraint(0,0.001);
3286 primVtxImprovedGamma+=v0GammaC;
3287 v0GammaC.SetProductionVertex(primVtxImprovedGamma);
3288
3289
3290 curGamma=curElecNeg+curElecPos;
3291 curGammaAt=curElecNegAt+curElecPosAt;
3292
3293 // invariant mass versus pt of K0short
3294
3295 Double_t chi2V0GammaC=100000.;
3296 if( v0GammaC.GetNDF() != 0) {
3297 chi2V0GammaC = v0GammaC.GetChi2()/v0GammaC.GetNDF();
3298 }else{
3299 cout<< "ERROR::v0K0C.GetNDF()" << endl;
3300 }
3301
3302 if(chi2V0GammaC<200 &&chi2V0GammaC>0 ){
3303 if(fHistograms != NULL){
3304 fHistograms->FillHistogram("ESD_RecalculateV0_InvMass",v0GammaC.GetMass());
3305 fHistograms->FillHistogram("ESD_RecalculateV0_Pt",v0GammaC.GetPt());
3306 fHistograms->FillHistogram("ESD_RecalculateV0_E_dEdxP",curElecNegAt.P(),negTrack->GetTPCsignal());
3307 fHistograms->FillHistogram("ESD_RecalculateV0_P_dEdxP",curElecPosAt.P(),posTrack->GetTPCsignal());
3308 fHistograms->FillHistogram("ESD_RecalculateV0_cpa",cpa);
3309 fHistograms->FillHistogram("ESD_RecalculateV0_dca",dca);
9c1cb6f7 3310 fHistograms->FillHistogram("ESD_RecalculateV0_normdistP",normdistP);
3311 fHistograms->FillHistogram("ESD_RecalculateV0_normdistN",normdistN);
3312
3313 new((*fKFRecalculatedGammasTClone)[fKFRecalculatedGammasTClone->GetEntriesFast()]) AliKFParticle(v0GammaC);
3314 fElectronRecalculatedv1.push_back(iTracks);
3315 fElectronRecalculatedv2.push_back(jTracks);
3316 }
3317 }
3318 }
3319 }
3320
3321 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFRecalculatedGammasTClone->GetEntriesFast();firstGammaIndex++){
3322 for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFRecalculatedGammasTClone->GetEntriesFast();secondGammaIndex++){
3323 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFRecalculatedGammasTClone->At(firstGammaIndex);
3324 AliKFParticle * twoGammaDecayCandidateDaughter1 = (AliKFParticle *)fKFRecalculatedGammasTClone->At(secondGammaIndex);
3325
3326 if(fElectronRecalculatedv1[firstGammaIndex]==fElectronRecalculatedv1[secondGammaIndex]){
3327 continue;
3328 }
3329 if( fElectronRecalculatedv2[firstGammaIndex]==fElectronRecalculatedv2[secondGammaIndex]){
3330 continue;
3331 }
3332
3333 AliKFParticle twoGammaCandidate(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
3334 fHistograms->FillHistogram("ESD_RecalculateGG_InvMass",twoGammaCandidate.GetMass());
3335 fHistograms->FillHistogram("ESD_RecalculateGG_InvMass_vs_Pt",twoGammaCandidate.GetMass(),twoGammaCandidate.GetPt());
3336 }
3337 }
3338}
01b7fdcc 3339void AliAnalysisTaskGammaConversion::CalculateJetCone(Int_t gammaIndex){
6c84d371 3340 // CaculateJetCone
a0b94e5c 3341
d7d7e825 3342 Double_t cone;
3343 Double_t coneSize=0.3;
3344 Double_t ptJet=0;
a0b94e5c 3345
6c84d371 3346 // AliKFParticle * currentGamma = &fKFReconstructedGammas[gammaIndex];
3347 AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
98778c17 3348
01b7fdcc 3349 TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
a0b94e5c 3350
6c84d371 3351 AliESDtrack* leadingCharged = (AliESDtrack*)(fChargedParticles->At(fLeadingChargedIndex));
98778c17 3352
d7d7e825 3353 Double_t momLeadingCharged[3];
3354 leadingCharged->GetConstrainedPxPyPz(momLeadingCharged);
a0b94e5c 3355
01b7fdcc 3356 TVector3 momentumVectorLeadingCharged(momLeadingCharged[0],momLeadingCharged[1],momLeadingCharged[2]);
a0b94e5c 3357
01b7fdcc 3358 Double_t phi1=momentumVectorLeadingCharged.Phi();
3359 Double_t eta1=momentumVectorLeadingCharged.Eta();
3360 Double_t phi3=momentumVectorCurrentGamma.Phi();
a0b94e5c 3361
6c84d371 3362 for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
3363 AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
d7d7e825 3364 Int_t chId = fChargedParticlesId[iCh];
3365 if(fLeadingChargedIndex==chId || fLeadingChargedIndex==chId) continue;
3366 Double_t mom[3];
3367 curTrack->GetConstrainedPxPyPz(mom);
01b7fdcc 3368 TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
3369 Double_t phi2=momentumVectorChargedParticle.Phi();
3370 Double_t eta2=momentumVectorChargedParticle.Eta();
a0b94e5c 3371
3372
d7d7e825 3373 cone=100.;
3374 if( TMath::Abs(phi2 - phi1) <= ( TMath::TwoPi()-coneSize) ){
3375 cone = TMath::Sqrt( TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-phi1),2) );
3376 }else{
3377 if( (phi2 - phi1)> TMath::TwoPi()-coneSize ){
3378 cone = TMath::Sqrt( TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-TMath::TwoPi()-phi1),2) );
3379 }
3380 if( (phi2 - phi1)< -(TMath::TwoPi()-coneSize) ){
3381 cone = TMath::Sqrt( TMath::Power((eta2-eta1),2)+ TMath::Power((phi2+TMath::TwoPi()-phi1),2) );
3382 }
3383 }
a0b94e5c 3384
01b7fdcc 3385 if(cone <coneSize&& momentumVectorChargedParticle.Pt()>fMinPtJetCone ){
3386 ptJet+= momentumVectorChargedParticle.Pt();
3387 Double_t ffzHdrGam = momentumVectorChargedParticle.Pt()/momentumVectorCurrentGamma.Pt();
3388 Double_t imbalanceHdrGam=-momentumVectorChargedParticle.Dot(momentumVectorCurrentGamma)/momentumVectorCurrentGamma.Mag2();
d7d7e825 3389 fHistograms->FillHistogram("ESD_FFzHdrGam",ffzHdrGam);
3390 fHistograms->FillHistogram("ESD_ImbalanceHdrGam",imbalanceHdrGam);
a0b94e5c 3391
d7d7e825 3392 }
a0b94e5c 3393
d7d7e825 3394 Double_t dphiHdrGam=phi3-phi2;
3395 if ( dphiHdrGam < (-TMath::PiOver2())){
3396 dphiHdrGam+=(TMath::TwoPi());
3397 }
a0b94e5c 3398
d7d7e825 3399 if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
3400 dphiHdrGam-=(TMath::TwoPi());
3401 }
a0b94e5c 3402
01b7fdcc 3403 if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
d7d7e825 3404 fHistograms->FillHistogram("ESD_dphiHdrGamIsolated",dphiHdrGam);
3405 }
3406 }//track loop
a0b94e5c 3407
3408
d7d7e825 3409}
3410
3411Double_t AliAnalysisTaskGammaConversion::GetMinimumDistanceToCharge(Int_t indexHighestPtGamma){
3412 // GetMinimumDistanceToCharge
a0b94e5c 3413
d7d7e825 3414 Double_t fIsoMin=100.;
3415 Double_t ptLeadingCharged=-1.;
98778c17 3416
3417 fLeadingChargedIndex=-1;
a0b94e5c 3418
6c84d371 3419 AliKFParticle * gammaHighestPt = (AliKFParticle*)fKFReconstructedGammasTClone->At(indexHighestPtGamma);
01b7fdcc 3420 TVector3 momentumVectorgammaHighestPt(gammaHighestPt->GetPx(),gammaHighestPt->GetPy(),gammaHighestPt->GetPz());
a0b94e5c 3421
01b7fdcc 3422 Double_t phi1=momentumVectorgammaHighestPt.Phi();
3423 Double_t eta1=momentumVectorgammaHighestPt.Eta();
a0b94e5c 3424
6c84d371 3425 for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
3426 AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
d7d7e825 3427 Int_t chId = fChargedParticlesId[iCh];
3428 if(fElectronv1[indexHighestPtGamma]==chId || fElectronv2[indexHighestPtGamma]==chId) continue;
3429 Double_t mom[3];
3430 curTrack->GetConstrainedPxPyPz(mom);
01b7fdcc 3431 TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
3432 Double_t phi2=momentumVectorChargedParticle.Phi();
3433 Double_t eta2=momentumVectorChargedParticle.Eta();
d7d7e825 3434 Double_t iso=pow( (pow( (eta1-eta2),2)+ pow((phi1-phi2),2)),0.5 );
a0b94e5c 3435
01b7fdcc 3436 if(momentumVectorChargedParticle.Pt()>fMinPtIsoCone ){
d7d7e825 3437 if (iso<fIsoMin){
3438 fIsoMin=iso;
3439 }
3440 }
a0b94e5c 3441
d7d7e825 3442 Double_t dphiHdrGam=phi1-phi2;
3443 if ( dphiHdrGam < (-TMath::PiOver2())){
3444 dphiHdrGam+=(TMath::TwoPi());
3445 }
a0b94e5c 3446
d7d7e825 3447 if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
3448 dphiHdrGam-=(TMath::TwoPi());
3449 }
01b7fdcc 3450 if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
d7d7e825 3451 fHistograms->FillHistogram("ESD_dphiHdrGam",dphiHdrGam);
3452 }
a0b94e5c 3453
d7d7e825 3454 if (dphiHdrGam>0.9*TMath::Pi() && dphiHdrGam<1.1*TMath::Pi()){
a0b94e5c 3455 if (momentumVectorChargedParticle.Pt()> ptLeadingCharged && momentumVectorChargedParticle.Pt()>0.1*momentumVectorgammaHighestPt.Pt()){
01b7fdcc 3456 ptLeadingCharged=momentumVectorChargedParticle.Pt();
d7d7e825 3457 fLeadingChargedIndex=iCh;
3458 }
3459 }
a0b94e5c 3460
d7d7e825 3461 }//track loop
3462 fHistograms->FillHistogram("ESD_MinimumIsoDistance",fIsoMin);
3463 return fIsoMin;
a0b94e5c 3464
d7d7e825 3465}
3466
a0b94e5c 3467Int_t AliAnalysisTaskGammaConversion::GetIndexHighestPtGamma(){
3468 //GetIndexHighestPtGamma
3469
d7d7e825 3470 Int_t indexHighestPtGamma=-1;
01b7fdcc 3471 //Double_t
3472 fGammaPtHighest = -100.;
a0b94e5c 3473
6c84d371 3474 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
3475 AliKFParticle * gammaHighestPtCandidate = (AliKFParticle*)fKFReconstructedGammasTClone->At(firstGammaIndex);
a0b94e5c 3476 TVector3 momentumVectorgammaHighestPtCandidate(gammaHighestPtCandidate->GetPx(),gammaHighestPtCandidate->GetPy(),gammaHighestPtCandidate->GetPz());
3477 if (momentumVectorgammaHighestPtCandidate.Pt() > fGammaPtHighest){
3478 fGammaPtHighest=momentumVectorgammaHighestPtCandidate.Pt();
3479 //gammaHighestPt = gammaHighestPtCandidate;
3480 indexHighestPtGamma=firstGammaIndex;
3481 }
d7d7e825 3482 }
a0b94e5c 3483
d7d7e825 3484 return indexHighestPtGamma;
a0b94e5c 3485
d7d7e825 3486}
3487
3488
3489void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)
3490{
3491 // Terminate analysis
3492 //
3493 AliDebug(1,"Do nothing in Terminate");
3494}
3495
3496void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()
3497{
332f1f44 3498
63477d82 3499 if(fKFCreateAOD) {
3500
3501 //AOD
3502 if(!fAODGamma) fAODGamma = new TClonesArray(fOutputAODClassName.Data(), 0);
3503 else fAODGamma->Delete();
3504 fAODGamma->SetName(Form("%s_gamma", fAODBranchName.Data()));
3505
3506 if(!fAODPi0) fAODPi0 = new TClonesArray(fOutputAODClassName.Data(), 0);
3507 else fAODPi0->Delete();
3508 fAODPi0->SetName(Form("%s_Pi0", fAODBranchName.Data()));
3509
3510 if(!fAODOmega) fAODOmega = new TClonesArray(fOutputAODClassName.Data(), 0);
3511 else fAODOmega->Delete();
3512 fAODOmega->SetName(Form("%s_Omega", fAODBranchName.Data()));
3513
3514 //If delta AOD file name set, add in separate file. Else add in standard aod file.
3515 if(GetDeltaAODFileName().Length() > 0) {
3516 AddAODBranch("TClonesArray", &fAODGamma, GetDeltaAODFileName().Data());
3517 AddAODBranch("TClonesArray", &fAODPi0, GetDeltaAODFileName().Data());
3518 AddAODBranch("TClonesArray", &fAODOmega, GetDeltaAODFileName().Data());
3519 AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(GetDeltaAODFileName().Data());
3520 } else {
3521 AddAODBranch("TClonesArray", &fAODGamma);
3522 AddAODBranch("TClonesArray", &fAODPi0);
3523 AddAODBranch("TClonesArray", &fAODOmega);
3524 }
332f1f44 3525 }
04bf4381 3526
d7d7e825 3527 // Create the output container
3528 if(fOutputContainer != NULL){
3529 delete fOutputContainer;
3530 fOutputContainer = NULL;
3531 }
3532 if(fOutputContainer == NULL){
3533 fOutputContainer = new TList();
b58d3c74 3534 fOutputContainer->SetOwner(kTRUE);
d7d7e825 3535 }
3536
3537 //Adding the histograms to the output container
3538 fHistograms->GetOutputContainer(fOutputContainer);
3539
3540
3541 if(fWriteNtuple){
3542 if(fGammaNtuple == NULL){
3543 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);
3544 }
3545 if(fNeutralMesonNtuple == NULL){
3546 fNeutralMesonNtuple = new TNtuple("NeutralMesonNtuple","NeutralMesonNtuple","test");
3547 }
3548 TList * ntupleTList = new TList();
b58d3c74 3549 ntupleTList->SetOwner(kTRUE);
d7d7e825 3550 ntupleTList->SetName("Ntuple");
3551 ntupleTList->Add((TNtuple*)fGammaNtuple);
3552 fOutputContainer->Add(ntupleTList);
3553 }
3554
3555 fOutputContainer->SetName(GetName());
3556}
3557
77ac6f3e 3558Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(const TParticle* const daughter0, const TParticle* const daughter1) const{
d7d7e825 3559 //helper function
3560 TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());
3561 TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());
3562 return v3D0.Angle(v3D1);
3563}
3564
3565void AliAnalysisTaskGammaConversion::CheckV0Efficiency(){
3566 // see header file for documentation
5e55d806 3567
d7d7e825 3568 vector<Int_t> indexOfGammaParticle;
a0b94e5c 3569
d7d7e825 3570 fStack = fV0Reader->GetMCStack();
a0b94e5c 3571
d7d7e825 3572 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
3573 return; // aborts if the primary vertex does not have contributors.
3574 }
a0b94e5c 3575
d7d7e825 3576 for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {
3577 TParticle* particle = (TParticle *)fStack->Particle(iTracks);
3578 if(particle->GetPdgCode()==22){ //Gamma
3579 if(particle->GetNDaughters() >= 2){
3580 TParticle* electron=NULL;
3581 TParticle* positron=NULL;
3582 for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
3583 TParticle *tmpDaughter = fStack->Particle(daughterIndex);
3584 if(tmpDaughter->GetUniqueID() == 5){
3585 if(tmpDaughter->GetPdgCode() == 11){
3586 electron = tmpDaughter;
3587 }
3588 else if(tmpDaughter->GetPdgCode() == -11){
3589 positron = tmpDaughter;
3590 }
3591 }
3592 }
3593 if(electron!=NULL && positron!=0){
3594 if(electron->R()<160){
3595 indexOfGammaParticle.push_back(iTracks);
3596 }
3597 }
3598 }
3599 }
3600 }
a0b94e5c 3601
d7d7e825 3602 Int_t nFoundGammas=0;
3603 Int_t nNotFoundGammas=0;
a0b94e5c 3604
d7d7e825 3605 Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
3606 for(Int_t i=0;i<numberOfV0s;i++){
3607 fV0Reader->GetV0(i);
a0b94e5c 3608
d7d7e825 3609 if(fV0Reader->HasSameMCMother() == kFALSE){
3610 continue;
3611 }
a0b94e5c 3612
d7d7e825 3613 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
3614 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
a0b94e5c 3615
d7d7e825 3616 if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
3617 continue;
3618 }
3619 if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
3620 continue;
3621 }
a0b94e5c 3622
d7d7e825 3623 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
3624 //TParticle * v0Gamma = fV0Reader->GetMotherMCParticle();
3625 for(UInt_t mcIndex=0;mcIndex<indexOfGammaParticle.size();mcIndex++){
3626 if(negativeMC->GetFirstMother()==indexOfGammaParticle[mcIndex]){
3627 nFoundGammas++;
3628 }
3629 else{
3630 nNotFoundGammas++;
3631 }
3632 }
3633 }
3634 }
d7d7e825 3635}
3636
3637
3638void AliAnalysisTaskGammaConversion::ProcessGammaElectronsForChicAnalysis(){
3639 // see header file for documantation
a0b94e5c 3640
d7d7e825 3641 fESDEvent = fV0Reader->GetESDEvent();
a0b94e5c 3642
3643
3644 TClonesArray * vESDeNegTemp = new TClonesArray("AliESDtrack",0);
6c84d371 3645 TClonesArray * vESDePosTemp = new TClonesArray("AliESDtrack",0);
3646 TClonesArray * vESDxNegTemp = new TClonesArray("AliESDtrack",0);
3647 TClonesArray * vESDxPosTemp = new TClonesArray("AliESDtrack",0);
3648 TClonesArray * vESDeNegNoJPsi = new TClonesArray("AliESDtrack",0);
3649 TClonesArray * vESDePosNoJPsi = new TClonesArray("AliESDtrack",0);
a0b94e5c 3650
6c84d371 3651 /*
3652 vector <AliESDtrack*> vESDeNegTemp(0);
3653 vector <AliESDtrack*> vESDePosTemp(0);
3654 vector <AliESDtrack*> vESDxNegTemp(0);
3655 vector <AliESDtrack*> vESDxPosTemp(0);
3656 vector <AliESDtrack*> vESDeNegNoJPsi(0);
3657 vector <AliESDtrack*> vESDePosNoJPsi(0);
3658 */
a0b94e5c 3659
3660
d7d7e825 3661 fHistograms->FillTable("Table_Electrons",0);//Count number of Events
a0b94e5c 3662
d7d7e825 3663 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
3664 AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
a0b94e5c 3665
d7d7e825 3666 if(!curTrack){
3667 //print warning here
3668 continue;
3669 }
a0b94e5c 3670
d7d7e825 3671 double p[3];if(!curTrack->GetConstrainedPxPyPz(p))continue;
3672 double r[3];curTrack->GetConstrainedXYZ(r);
a0b94e5c 3673
d7d7e825 3674 TVector3 rXYZ(r);
a0b94e5c 3675
d7d7e825 3676 fHistograms->FillTable("Table_Electrons",4);//Count number of ESD tracks
a0b94e5c 3677
d7d7e825 3678 Bool_t flagKink = kTRUE;
3679 Bool_t flagTPCrefit = kTRUE;
3680 Bool_t flagTRDrefit = kTRUE;
3681 Bool_t flagITSrefit = kTRUE;
3682 Bool_t flagTRDout = kTRUE;
3683 Bool_t flagVertex = kTRUE;
a0b94e5c 3684
3685
d7d7e825 3686 //Cuts ---------------------------------------------------------------
a0b94e5c 3687
d7d7e825 3688 if(curTrack->GetKinkIndex(0) > 0){
3689 fHistograms->FillHistogram("Table_Electrons",5);//Count kink
3690 flagKink = kFALSE;
3691 }
a0b94e5c 3692
d7d7e825 3693 ULong_t trkStatus = curTrack->GetStatus();
a0b94e5c 3694
d7d7e825 3695 ULong_t tpcRefit = (trkStatus & AliESDtrack::kTPCrefit);
a0b94e5c 3696
d7d7e825 3697 if(!tpcRefit){
3698 fHistograms->FillHistogram("Table_Electrons",9);//Count not TPCrefit
3699 flagTPCrefit = kFALSE;
3700 }
a0b94e5c 3701
d7d7e825 3702 ULong_t itsRefit = (trkStatus & AliESDtrack::kITSrefit);
3703 if(!itsRefit){
3704 fHistograms->FillHistogram("Table_Electrons",10);//Count not ITSrefit
3705 flagITSrefit = kFALSE;
3706 }
a0b94e5c 3707
d7d7e825 3708 ULong_t trdRefit = (trkStatus & AliESDtrack::kTRDrefit);
a0b94e5c 3709
d7d7e825 3710 if(!trdRefit){
3711 fHistograms->FillHistogram("Table_Electrons",8); //Count not TRDrefit
3712 flagTRDrefit = kFALSE;
3713 }
a0b94e5c 3714
d7d7e825 3715 ULong_t trdOut = (trkStatus & AliESDtrack::kTRDout);
a0b94e5c 3716
d7d7e825 3717 if(!trdOut) {
3718 fHistograms->FillHistogram("Table_Electrons",7); //Count not TRDout
3719 flagTRDout = kFALSE;
3720 }
a0b94e5c 3721
d7d7e825 3722 double nsigmaToVxt = GetSigmaToVertex(curTrack);
a0b94e5c 3723
d7d7e825 3724 if(nsigmaToVxt > 3){
3725 fHistograms->FillHistogram("Table_Electrons",6); //Count Tracks with number of sigmas > 3
3726 flagVertex = kFALSE;
3727 }
a0b94e5c 3728
d7d7e825 3729 if(! (flagKink && flagTPCrefit && flagITSrefit && flagTRDrefit && flagTRDout && flagVertex ) ) continue;
3730 fHistograms->FillHistogram("Table_Electrons",11);//Count Tracks passed Cuts
a0b94e5c 3731
3732
d7d7e825 3733 Stat_t pid, weight;
3734 GetPID(curTrack, pid, weight);
a0b94e5c 3735
d7d7e825 3736 if(pid!=0){
3737 fHistograms->FillHistogram("Table_Electrons",12); //Count Tracks with pid != 0
3738 }
a0b94e5c 3739
d7d7e825 3740 if(pid == 0){
3741 fHistograms->FillHistogram("Table_Electrons",13); //Count Tracks with pid != 0
3742 }
a0b94e5c 3743
3744
3745
3746
a0b94e5c 3747
3748
d7d7e825 3749 TLorentzVector curElec;
3750 curElec.SetXYZM(p[0],p[1],p[2],fElectronMass);
a0b94e5c 3751
3752
113d8432 3753 if(fDoMCTruth){
3754 Int_t labelMC = TMath::Abs(curTrack->GetLabel());
3755 TParticle* curParticle = fStack->Particle(labelMC);
3756 if(curTrack->GetSign() > 0){
3757 if( pid == 0){
3758 fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
3759 fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
3760 }
3761 else{
3762 fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
3763 fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
3764 }
3765 }
3766 }
a0b94e5c 3767
3768
d7d7e825 3769 if(curTrack->GetSign() > 0){
a0b94e5c 3770
6c84d371 3771 // vESDxPosTemp.push_back(curTrack);
3772 new((*vESDxPosTemp)[vESDxPosTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
a0b94e5c 3773
d7d7e825 3774 if( pid == 0){
a0b94e5c 3775
d7d7e825 3776 fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
3777 fHistograms->FillHistogram("ESD_ElectronPosPt",curElec.Pt());
113d8432 3778 // fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
d7d7e825 3779 fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
113d8432 3780 // fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
6c84d371 3781 // vESDePosTemp.push_back(curTrack);
3782 new((*vESDePosTemp)[vESDePosTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
a0b94e5c 3783
d7d7e825 3784 }
a0b94e5c 3785
d7d7e825 3786 }
3787 else {
5e55d806 3788
6c84d371 3789 new((*vESDxNegTemp)[vESDxNegTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
a0b94e5c 3790
d7d7e825 3791 if( pid == 0){
a0b94e5c 3792
d7d7e825 3793 fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
3794 fHistograms->FillHistogram("ESD_ElectronNegPt",curElec.Pt());
d7d7e825 3795 fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
6c84d371 3796 new((*vESDeNegTemp)[vESDeNegTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
a0b94e5c 3797
d7d7e825 3798 }
a0b94e5c 3799
d7d7e825 3800 }
a0b94e5c 3801
d7d7e825 3802 }
a0b94e5c 3803
3804
d7d7e825 3805 Bool_t ePosJPsi = kFALSE;
3806 Bool_t eNegJPsi = kFALSE;
3807 Bool_t ePosPi0 = kFALSE;
3808 Bool_t eNegPi0 = kFALSE;
3809
3810 UInt_t iePosJPsi=0,ieNegJPsi=0,iePosPi0=0,ieNegPi0=0;
a0b94e5c 3811
6c84d371 3812 for(Int_t iNeg=0; iNeg < vESDeNegTemp->GetEntriesFast(); iNeg++){
3813 if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetPdgCode() == 11)
3814 if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetMother(0) > -1){
3815 Int_t labelMother = fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetMother(0);
d7d7e825 3816 TParticle* partMother = fStack ->Particle(labelMother);
3817 if (partMother->GetPdgCode() == 111){
3818 ieNegPi0 = iNeg;
3819 eNegPi0 = kTRUE;
3820 }
3821 if(partMother->GetPdgCode() == 443){ //Mother JPsi
3822 fHistograms->FillTable("Table_Electrons",14);
3823 ieNegJPsi = iNeg;
3824 eNegJPsi = kTRUE;
3825 }
3826 else{
6c84d371 3827 // vESDeNegNoJPsi.push_back(vESDeNegTemp[iNeg]);
3828 new((*vESDeNegNoJPsi)[vESDeNegNoJPsi->GetEntriesFast()]) AliESDtrack(*(AliESDtrack*)(vESDeNegTemp->At(iNeg)));
d7d7e825 3829 // cout<<"ESD No Positivo JPsi "<<endl;
3830 }
a0b94e5c 3831
d7d7e825 3832 }
3833 }
a0b94e5c 3834
6c84d371 3835 for(Int_t iPos=0; iPos < vESDePosTemp->GetEntriesFast(); iPos++){
3836 if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetPdgCode() == -11)
3837 if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetMother(0) > -1){
3838 Int_t labelMother = fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetMother(0);
d7d7e825 3839 TParticle* partMother = fStack ->Particle(labelMother);
3840 if (partMother->GetPdgCode() == 111){
3841 iePosPi0 = iPos;
3842 ePosPi0 = kTRUE;
3843 }
3844 if(partMother->GetPdgCode() == 443){ //Mother JPsi
3845 fHistograms->FillTable("Table_Electrons",15);
3846 iePosJPsi = iPos;
3847 ePosJPsi = kTRUE;
3848 }
3849 else{
6c84d371 3850 // vESDePosNoJPsi.push_back(vESDePosTemp[iPos]);
3851 new((*vESDePosNoJPsi)[vESDePosNoJPsi->GetEntriesFast()]) AliESDtrack(*(AliESDtrack*)(vESDePosTemp->At(iPos)));
d7d7e825 3852 // cout<<"ESD No Negativo JPsi "<<endl;
3853 }
a0b94e5c 3854
d7d7e825 3855 }
3856 }
3857
3858 if( eNegJPsi && ePosJPsi ){
3859 TVector3 tempeNegV,tempePosV;
6c84d371 3860 tempeNegV.SetXYZ(((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Px(),((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Py(),((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Pz());
3861 tempePosV.SetXYZ(((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Px(),((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Py(),((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Pz());
d7d7e825 3862 fHistograms->FillTable("Table_Electrons",16);
3863 fHistograms->FillHistogram("ESD_ElectronPosNegJPsiAngle",tempeNegV.Angle(tempePosV));
6c84d371 3864 fHistograms->FillHistogram("MC_ElectronPosNegJPsiAngle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->GetLabel())),
3865 fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->GetLabel()))));
d7d7e825 3866 }
3867
3868 if( eNegPi0 && ePosPi0 ){
3869 TVector3 tempeNegV,tempePosV;
6c84d371 3870 tempeNegV.SetXYZ(((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Px(),((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Py(),((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Pz());
3871 tempePosV.SetXYZ(((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Px(),((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Py(),((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Pz());
d7d7e825 3872 fHistograms->FillHistogram("ESD_ElectronPosNegPi0Angle",tempeNegV.Angle(tempePosV));
6c84d371 3873 fHistograms->FillHistogram("MC_ElectronPosNegPi0Angle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->GetLabel())),
3874 fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->GetLabel()))));
d7d7e825 3875 }
a0b94e5c 3876
3877
d7d7e825 3878 FillAngle("ESD_eNegePosAngleBeforeCut",GetTLorentzVector(vESDeNegTemp),GetTLorentzVector(vESDePosTemp));
a0b94e5c 3879
6c84d371 3880 CleanWithAngleCuts(*vESDeNegTemp,*vESDePosTemp,*fKFReconstructedGammasTClone);
d7d7e825 3881
6c84d371 3882 // vector <TLorentzVector> vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectron);
3883 // vector <TLorentzVector> vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectron);
a0b94e5c 3884
6c84d371 3885 TClonesArray vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectronTClone);
3886 TClonesArray vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectronTClone);
a0b94e5c 3887
3888
d7d7e825 3889 FillAngle("ESD_eNegePosAngleAfterCut",vCurrentTLVeNeg,vCurrentTLVePos);
a0b94e5c 3890
3891
3892
3893
d7d7e825 3894 //FillAngle("ESD_eNegePosAngleAfterCut",CurrentTLVeNeg,CurrentTLVePos);
a0b94e5c 3895
3896
d7d7e825 3897 FillElectronInvMass("ESD_InvMass_ePluseMinus",vCurrentTLVeNeg,vCurrentTLVePos);
3898 FillElectronInvMass("ESD_InvMass_xPlusxMinus",GetTLorentzVector(vESDxNegTemp),GetTLorentzVector(vESDxPosTemp));
a0b94e5c 3899
3900
3901
6c84d371 3902 FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusChiC","ESD_InvMass_GammaePluseMinusChiCDiff",*fKFReconstructedGammasCutTClone,vCurrentTLVeNeg,vCurrentTLVePos);
a0b94e5c 3903
d7d7e825 3904 FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusPi0","ESD_InvMass_GammaePluseMinusPi0Diff",
6c84d371 3905 *fKFReconstructedGammasCutTClone,vCurrentTLVeNeg,vCurrentTLVePos);
a0b94e5c 3906
d7d7e825 3907 //BackGround
a0b94e5c 3908
d7d7e825 3909 //Like Sign e+e-
3910 ElectronBackground("ESD_ENegBackground",vCurrentTLVeNeg);
3911 ElectronBackground("ESD_EPosBackground",vCurrentTLVePos);
3912 ElectronBackground("ESD_EPosENegBackground",vCurrentTLVeNeg);
3913 ElectronBackground("ESD_EPosENegBackground",vCurrentTLVePos);
a0b94e5c 3914
d7d7e825 3915 // Like Sign e+e- no JPsi
3916 ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDeNegNoJPsi));
3917 ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDePosNoJPsi));
a0b94e5c 3918
d7d7e825 3919 //Mixed Event
a0b94e5c 3920
6c84d371 3921 if( fCurrentEventPosElectronTClone->GetEntriesFast() > 0 && fCurrentEventNegElectronTClone->GetEntriesFast() > 0 && fKFReconstructedGammasCutTClone->GetEntriesFast() > 0 ){
d7d7e825 3922 FillGammaElectronInvMass("ESD_EPosENegGammaBackgroundMX","ESD_EPosENegGammaBackgroundMXDiff",
6c84d371 3923 *fKFReconstructedGammasCutTClone,*fPreviousEventTLVNegElectronTClone,*fPreviousEventTLVPosElectronTClone);
3924 *fPreviousEventTLVNegElectronTClone = vCurrentTLVeNeg;
3925 *fPreviousEventTLVPosElectronTClone = vCurrentTLVePos;
a0b94e5c 3926
d7d7e825 3927 }
a0b94e5c 3928
d7d7e825 3929 /*
3930 //Photons P
3931 Double_t vtx[3];
3932 vtx[0]=0;vtx[1]=0;vtx[2]=0;
3933 for(UInt_t i=0;i<fKFReconstructedGammasChic.size();i++){
a0b94e5c 3934
d7d7e825 3935 // if(fMCGammaChicTempCut[i]->GetMother(0) < 0) continue;
a0b94e5c 3936
3937
3938
d7d7e825 3939 Int_t tempLabel = fStack->Particle(fMCGammaChicTempCut[i]->GetMother(0))->GetPdgCode();
3940 // cout<<" Label Pedro Gonzalez " <<tempLabel <<endl;
a0b94e5c 3941
d7d7e825 3942 // cout<<" Label Distance"<<fKFReconstructedGammasChic[i].GetDistanceFromVertex(vtx)<<endl;
a0b94e5c 3943
d7d7e825 3944 if( tempLabel == 10441 || tempLabel == 20443 || tempLabel == 445 )
a0b94e5c 3945
d7d7e825 3946 fHistograms->FillHistogram("ESD_PhotonsMomentum",fKFReconstructedGammasChic[i].GetMomentum());
a0b94e5c 3947
3948
d7d7e825 3949 }
a0b94e5c 3950
3951
d7d7e825 3952 */
1e7846f4 3953
3954
3955 vESDeNegTemp->Delete();
3956 vESDePosTemp->Delete();
3957 vESDxNegTemp->Delete();
3958 vESDxPosTemp->Delete();
3959 vESDeNegNoJPsi->Delete();
3960 vESDePosNoJPsi->Delete();
3961
3962 delete vESDeNegTemp;
3963 delete vESDePosTemp;
3964 delete vESDxNegTemp;
3965 delete vESDxPosTemp;
3966 delete vESDeNegNoJPsi;
3967 delete vESDePosNoJPsi;
d7d7e825 3968}
3969
6c84d371 3970/*
3971 void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,vector <TLorentzVector> tlVeNeg, vector <TLorentzVector> tlVePos){
d7d7e825 3972 //see header file for documentation
3973 for( UInt_t iNeg=0; iNeg < tlVeNeg.size(); iNeg++){
6c84d371 3974 for (UInt_t iPos=0; iPos < tlVePos.size(); iPos++){
3975 fHistograms->FillHistogram(histoName.Data(),tlVeNeg[iNeg].Vect().Angle(tlVePos[iPos].Vect()));
3976 }
3977 }
3978 }
3979*/
3980void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,TClonesArray const tlVeNeg, TClonesArray const tlVePos){
3981 //see header file for documentation
3982 for( Int_t iNeg=0; iNeg < tlVeNeg.GetEntriesFast(); iNeg++){
3983 for (Int_t iPos=0; iPos < tlVePos.GetEntriesFast(); iPos++){
3984 fHistograms->FillHistogram(histoName.Data(),((TLorentzVector*)(tlVeNeg.At(iNeg)))->Vect().Angle(((TLorentzVector*)(tlVePos.At(iPos)))->Vect()));
d7d7e825 3985 }
3986 }
3987}
6c84d371 3988void AliAnalysisTaskGammaConversion::FillElectronInvMass(TString histoName, TClonesArray const eNeg, TClonesArray const ePos){
d7d7e825 3989 //see header file for documentation
6c84d371 3990 for( Int_t n=0; n < eNeg.GetEntriesFast(); n++){
3991 TLorentzVector en = (*(TLorentzVector*)(eNeg.At(n)));
3992 for (Int_t p=0; p < ePos.GetEntriesFast(); p++){
3993 TLorentzVector ep = (*(TLorentzVector*)(ePos.At(p)));
d7d7e825 3994 TLorentzVector np = ep + en;
3995 fHistograms->FillHistogram(histoName.Data(),np.M());
3996 }
3997 }
d7d7e825 3998}
3999
6c84d371 4000void AliAnalysisTaskGammaConversion::FillGammaElectronInvMass(TString histoMass,TString histoDiff,TClonesArray const fKFGammas,
4001 TClonesArray const tlVeNeg,TClonesArray const tlVePos)
d7d7e825 4002{
4003 //see header file for documentation
a0b94e5c 4004
6c84d371 4005 for( Int_t iNeg=0; iNeg < tlVeNeg.GetEntriesFast(); iNeg++ ){
a0b94e5c 4006
6c84d371 4007 for (Int_t iPos=0; iPos < tlVePos.GetEntriesFast(); iPos++){
a0b94e5c 4008
6c84d371 4009 TLorentzVector xy = *((TLorentzVector *)(tlVePos.At(iPos))) + *((TLorentzVector *)(tlVeNeg.At(iNeg)));
a0b94e5c 4010
6c84d371 4011 for (Int_t iGam=0; iGam < fKFGammas.GetEntriesFast(); iGam++){
a0b94e5c 4012
6c84d371 4013 // AliKFParticle * gammaCandidate = &fKFGammas[iGam];
4014 AliKFParticle * gammaCandidate = (AliKFParticle *)(fKFGammas.At(iGam));
d7d7e825 4015 TLorentzVector g;
a0b94e5c 4016
d7d7e825 4017 g.SetXYZM(gammaCandidate->GetPx(),gammaCandidate->GetPy(),gammaCandidate->GetPz(),fGammaMass);
4018 TLorentzVector xyg = xy + g;
4019 fHistograms->FillHistogram(histoMass.Data(),xyg.M());
4020 fHistograms->FillHistogram(histoDiff.Data(),(xyg.M()-xy.M()));
4021 }
4022 }
4023 }
a0b94e5c 4024
d7d7e825 4025}
6c84d371 4026void AliAnalysisTaskGammaConversion::ElectronBackground(TString hBg, TClonesArray e)
d7d7e825 4027{
4028 // see header file for documentation
6c84d371 4029 for(Int_t i=0; i < e.GetEntriesFast(); i++)
d7d7e825 4030 {
6c84d371 4031 for (Int_t j=i+1; j < e.GetEntriesFast(); j++)
d7d7e825 4032 {
6c84d371 4033 TLorentzVector ee = (*(TLorentzVector*)(e.At(i))) + (*(TLorentzVector*)(e.At(j)));
a0b94e5c 4034
d7d7e825 4035 fHistograms->FillHistogram(hBg.Data(),ee.M());
4036 }
4037 }
4038}
4039
4040
6c84d371 4041void AliAnalysisTaskGammaConversion::CleanWithAngleCuts(TClonesArray const negativeElectrons,
4042 TClonesArray const positiveElectrons,
4043 TClonesArray const gammas){
d7d7e825 4044 // see header file for documentation
a0b94e5c 4045
6c84d371 4046 UInt_t sizeN = negativeElectrons.GetEntriesFast();
4047 UInt_t sizeP = positiveElectrons.GetEntriesFast();
4048 UInt_t sizeG = gammas.GetEntriesFast();
a0b94e5c 4049
4050
4051
d7d7e825 4052 vector <Bool_t> xNegBand(sizeN);
4053 vector <Bool_t> xPosBand(sizeP);
4054 vector <Bool_t> gammaBand(sizeG);
a0b94e5c 4055
4056
d7d7e825 4057 for(UInt_t iNeg=0; iNeg < sizeN; iNeg++) xNegBand[iNeg]=kTRUE;
4058 for(UInt_t iPos=0; iPos < sizeP; iPos++) xPosBand[iPos]=kTRUE;
4059 for(UInt_t iGam=0; iGam < sizeG; iGam++) gammaBand[iGam]=kTRUE;
d7d7e825 4060
a0b94e5c 4061
4062 for(UInt_t iPos=0; iPos < sizeP; iPos++){
4063
6c84d371 4064 Double_t aP[3];
4065 ((AliESDtrack*)(positiveElectrons.At(iPos)))->GetConstrainedPxPyPz(aP);
a0b94e5c 4066
d7d7e825 4067 TVector3 ePosV(aP[0],aP[1],aP[2]);
a0b94e5c 4068
d7d7e825 4069 for(UInt_t iNeg=0; iNeg < sizeN; iNeg++){
a0b94e5c 4070
6c84d371 4071 Double_t aN[3];
4072 ((AliESDtrack*)(negativeElectrons.At(iNeg)))->GetConstrainedPxPyPz(aN);
d7d7e825 4073 TVector3 eNegV(aN[0],aN[1],aN[2]);
a0b94e5c 4074
d7d7e825 4075 if(ePosV.Angle(eNegV) < 0.05){ //e+e- from gamma
4076 xPosBand[iPos]=kFALSE;
4077 xNegBand[iNeg]=kFALSE;
4078 }
a0b94e5c 4079
d7d7e825 4080 for(UInt_t iGam=0; iGam < sizeG; iGam++){
6c84d371 4081 AliKFParticle* gammaCandidate = (AliKFParticle*)gammas.At(iGam);
d7d7e825 4082 TVector3 gammaCandidateVector(gammaCandidate->Px(),gammaCandidate->Py(),gammaCandidate->Pz());
4083 if(ePosV.Angle(gammaCandidateVector) < 0.05 || eNegV.Angle(gammaCandidateVector) < 0.05)
4084 gammaBand[iGam]=kFALSE;
4085 }
4086 }
4087 }
a0b94e5c 4088
4089
4090
4091
d7d7e825 4092 for(UInt_t iPos=0; iPos < sizeP; iPos++){
4093 if(xPosBand[iPos]){
6c84d371 4094 new((*fCurrentEventPosElectronTClone)[fCurrentEventPosElectronTClone->GetEntriesFast()]) AliESDtrack((*(AliESDtrack*)(positiveElectrons.At(iPos))));
4095 // fCurrentEventPosElectron.push_back(positiveElectrons[iPos]);
d7d7e825 4096 }
4097 }
4098 for(UInt_t iNeg=0;iNeg < sizeN; iNeg++){
4099 if(xNegBand[iNeg]){
6c84d371 4100 new((*fCurrentEventNegElectronTClone)[fCurrentEventNegElectronTClone->GetEntriesFast()]) AliESDtrack((*(AliESDtrack*)(negativeElectrons.At(iNeg))));
4101 // fCurrentEventNegElectron.push_back(negativeElectrons[iNeg]);
d7d7e825 4102 }
4103 }
4104 for(UInt_t iGam=0; iGam < sizeG; iGam++){
4105 if(gammaBand[iGam]){
6c84d371 4106 new((*fKFReconstructedGammasCutTClone)[fKFReconstructedGammasCutTClone->GetEntriesFast()]) AliKFParticle((*(AliKFParticle*)(gammas.At(iGam))));
4107 //fKFReconstructedGammasCut.push_back(*(AliKFParticle*)gammas->At(iGam));
d7d7e825 4108 }
4109 }
4110}
4111
4112
77ac6f3e 4113void AliAnalysisTaskGammaConversion::GetPID(const AliESDtrack *track, Stat_t &pid, Stat_t &weight)
d7d7e825 4114{
4115 // see header file for documentation
4116 pid = -1;
4117 weight = -1;
a0b94e5c 4118
d7d7e825 4119 double wpart[5];
4120 double wpartbayes[5];
a0b94e5c 4121
d7d7e825 4122 //get probability of the diffenrent particle types
4123 track->GetESDpid(wpart);
a0b94e5c 4124
d7d7e825 4125 // Tentative particle type "concentrations"
4126 double c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
a0b94e5c 4127
d7d7e825 4128 //Bayes' formula
4129 double rcc = 0.;
4130 for (int i = 0; i < 5; i++)
4131 {
4132 rcc+=(c[i] * wpart[i]);
4133 }
a0b94e5c 4134
4135
4136
d7d7e825 4137 for (int i=0; i<5; i++) {
4a6157dc 4138 if( rcc>0 || rcc<0){//Kenneth: not sure if the rcc<0 should be there, this is from fixing a coding violation where we are not allowed to say: rcc!=0 (RC19)
d7d7e825 4139 wpartbayes[i] = c[i] * wpart[i] / rcc;
4140 }
4141 }
a0b94e5c 4142
4143
4144
d7d7e825 4145 Float_t max=0.;
4146 int ipid=-1;
4147 //find most probable particle in ESD pid
4148 //0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton
4149 for (int i = 0; i < 5; i++)
4150 {
4151 if (wpartbayes[i] > max)
4152 {
a0b94e5c 4153 ipid = i;
4154 max = wpartbayes[i];
d7d7e825 4155 }
4156 }
a0b94e5c 4157
d7d7e825 4158 pid = ipid;
4159 weight = max;
4160}
77ac6f3e 4161double AliAnalysisTaskGammaConversion::GetSigmaToVertex(const AliESDtrack* t)
d7d7e825 4162{
4163 // Calculates the number of sigma to the vertex.
a0b94e5c 4164
d7d7e825 4165 Float_t b[2];
4166 Float_t bRes[2];
4167 Float_t bCov[3];
4168 t->GetImpactParameters(b,bCov);
4169 if (bCov[0]<=0 || bCov[2]<=0) {
4170 AliDebug(1, "Estimated b resolution lower or equal zero!");
4171 bCov[0]=0; bCov[2]=0;
4172 }
4173 bRes[0] = TMath::Sqrt(bCov[0]);
4174 bRes[1] = TMath::Sqrt(bCov[2]);
a0b94e5c 4175
d7d7e825 4176 // -----------------------------------
4177 // How to get to a n-sigma cut?
4178 //
4179 // The accumulated statistics from 0 to d is
4180 //
4181 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
4182 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
4183 //
4184 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
4185 // Can this be expressed in a different way?
a0b94e5c 4186
d7d7e825 4187 if (bRes[0] == 0 || bRes[1] ==0)
4188 return -1;
a0b94e5c 4189
d7d7e825 4190 double d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
a0b94e5c 4191
d7d7e825 4192 // stupid rounding problem screws up everything:
4193 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
4194 if (TMath::Exp(-d * d / 2) < 1e-10)
4195 return 1000;
a0b94e5c 4196
4197
d7d7e825 4198 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
4199 return d;
4200}
6c84d371 4201
4202//vector <TLorentzVector> AliAnalysisTaskGammaConversion::GetTLorentzVector(vector <AliESDtrack*> esdTrack){
4203TClonesArray AliAnalysisTaskGammaConversion::GetTLorentzVector(TClonesArray *const esdTrack){
01b7fdcc 4204 //Return TLoresntz vector of track?
6c84d371 4205 // vector <TLorentzVector> tlVtrack(0);
4206 TClonesArray array("TLorentzVector",0);
a0b94e5c 4207
6c84d371 4208 for(Int_t itrack=0; itrack < esdTrack->GetEntriesFast(); itrack++){
4209 double p[3];
4210 //esdTrack[itrack]->GetConstrainedPxPyPz(p);
4211 ((AliESDtrack*)(esdTrack->At(itrack)))->GetConstrainedPxPyPz(p);
d7d7e825 4212 TLorentzVector currentTrack;
01b7fdcc 4213 currentTrack.SetXYZM(p[0],p[1],p[2],fElectronMass);
6c84d371 4214 new((array)[array.GetEntriesFast()]) TLorentzVector(currentTrack);
4215 // tlVtrack.push_back(currentTrack);
d7d7e825 4216 }
a0b94e5c 4217
6c84d371 4218 return array;
a0b94e5c 4219
6c84d371 4220 // return tlVtrack;
d7d7e825 4221}
77ac6f3e 4222Int_t AliAnalysisTaskGammaConversion::GetProcessType(const AliMCEvent * mcEvt) {
3c45d101 4223
4224 // Determine if the event was generated with pythia or phojet and return the process type
4225
4226 // Check if mcEvt is fine
4227 if (!mcEvt) {
4228 AliFatal("NULL mc event");
4229 }
4230
4231 // Determine if it was a pythia or phojet header, and return the correct process type
4232 AliGenPythiaEventHeader * headPy = 0;
4233 AliGenDPMjetEventHeader * headPho = 0;
4234 AliGenEventHeader * htmp = mcEvt->GenEventHeader();
4235 if(!htmp) {
4236 AliFatal("Cannot Get MC Header!!");
4237 }
4238 if( TString(htmp->IsA()->GetName()) == "AliGenPythiaEventHeader") {
4239 headPy = (AliGenPythiaEventHeader*) htmp;
4240 } else if (TString(htmp->IsA()->GetName()) == "AliGenDPMjetEventHeader") {
4241 headPho = (AliGenDPMjetEventHeader*) htmp;
4242 } else {
4243 AliError("Unknown header");
4244 }
4245
4246 // Determine process type
4247 if(headPy) {
4248 if(headPy->ProcessType() == 92 || headPy->ProcessType() == 93) {
4249 // single difractive
4250 return kProcSD;
4251 } else if (headPy->ProcessType() == 94) {
4252 // double diffractive
4253 return kProcDD;
4254 }
4255 else if(headPy->ProcessType() != 92 && headPy->ProcessType() != 93 && headPy->ProcessType() != 94) {
4256 // non difractive
4257 return kProcND;
4258 }
4259 } else if (headPho) {
4260 if(headPho->ProcessType() == 5 || headPho->ProcessType() == 6 ) {
4261 // single difractive
4262 return kProcSD;
4263 } else if (headPho->ProcessType() == 7) {
4264 // double diffractive
4265 return kProcDD;
4266 } else if(headPho->ProcessType() != 5 && headPho->ProcessType() != 6 && headPho->ProcessType() != 7 ) {
4267 // non difractive
4268 return kProcND;
4269 }
4270 }
4271
4272
4273 // no process type found?
4274 AliError(Form("Unknown header: %s", htmp->IsA()->GetName()));
4275 return kProcUnknown;
4276}
6746e1e1 4277
4278
4279Int_t AliAnalysisTaskGammaConversion::CalculateMultiplicityBin(){
4280 // Get Centrality bin
4281
4282 Int_t multiplicity = 0;
4283
4284 if ( fUseMultiplicity == 1 ) {
4285
4286 if (fMultiplicity>= 0 && fMultiplicity<= 5) multiplicity=1;
4287 if (fMultiplicity>= 6 && fMultiplicity<= 9) multiplicity=2;
4288 if (fMultiplicity>=10 && fMultiplicity<=14) multiplicity=3;
4289 if (fMultiplicity>=15 && fMultiplicity<=22) multiplicity=4;
4290 if (fMultiplicity>=23 ) multiplicity=5;
4291
4292 }
4293 return multiplicity;
4294}