AliPhysics  2b88e80 (2b88e80)
AliAnalysisTaskHeavyNeutralMesonToGG.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2020, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: Friederike Bock, Sandro Wenzel *
5  * Version 1.0 *
6  * *
7  * *
8  * Permission to use, copy, modify and distribute this software and its *
9  * documentation strictly for non-commercial purposes is hereby granted *
10  * without fee, provided that the above copyright notice appears in all *
11  * copies and that both the copyright notice and this permission notice *
12  * appear in the supporting documentation. The authors make no claims *
13  * about the suitability of this software for any purpose. It is *
14  * provided "as is" without express or implied warranty. *
15  **************************************************************************/
16 
18 //----------------------------------------------------------------
19 // Class used to do analysis on conversion photons + calo photons
20 //----------------------------------------------------------------
22 #include "TChain.h"
23 #include "TTree.h"
24 #include "TBranch.h"
25 #include "TFile.h"
26 #include "TH1F.h"
27 #include "TH2F.h"
28 #include "TH3F.h"
29 #include "THnSparse.h"
30 #include "TCanvas.h"
31 #include "TNtuple.h"
32 #include "AliAnalysisTask.h"
33 #include "AliAnalysisManager.h"
34 #include "AliESDEvent.h"
35 #include "AliESDInputHandler.h"
36 #include "AliMCEventHandler.h"
37 #include "AliMCEvent.h"
38 #include "AliMCParticle.h"
39 #include "AliCentrality.h"
40 #include "AliESDVZERO.h"
41 #include "AliESDpid.h"
42 #include "AliAnalysisTaskHeavyNeutralMesonToGG.h"
43 #include "AliVParticle.h"
44 #include "AliESDtrack.h"
45 #include "AliESDtrackCuts.h"
46 #include "AliKFVertex.h"
47 #include "AliGenCocktailEventHeader.h"
48 #include "AliConversionAODBGHandlerRP.h"
49 #include "AliAODMCParticle.h"
50 #include "AliAODMCHeader.h"
51 #include "AliEventplane.h"
52 #include "AliAnalysisTaskEMCALClusterizeFast.h"
53 #include "AliAODEvent.h"
54 #include "AliAODInputHandler.h"
55 #include "AliESDEvent.h"
56 #include "AliESDInputHandler.h"
57 #include "AliInputEventHandler.h"
58 #include "AliCaloTrackMatcher.h"
59 #include <vector>
60 #include <map>
61 
62 ClassImp(AliAnalysisTaskHeavyNeutralMesonToGG)
63 
64 //________________________________________________________________________
65 AliAnalysisTaskHeavyNeutralMesonToGG::AliAnalysisTaskHeavyNeutralMesonToGG(): AliAnalysisTaskSE(),
66  fV0Reader(NULL),
67  fV0ReaderName("V0ReaderV1"),
68  fCorrTaskSetting(""),
69  fBGHandler(NULL),
70  fBGHandlerRP(NULL),
71  fBGClusHandler(NULL),
72  fBGClusHandlerRP(NULL),
73  fInputEvent(NULL),
74  fMCEvent(NULL),
75  fCutFolder(NULL),
76  fESDList(NULL),
77  fBackList(NULL),
78  fMotherList(NULL),
79  fTrueList(NULL),
80  fMCList(NULL),
81  fOutputContainer(NULL),
82  fReaderGammas(NULL),
83  fGammaCandidates(NULL),
84  fClusterCandidates(NULL),
85  fEventCutArray(NULL),
86  fEventCuts(NULL),
87  fCutArray(NULL),
88  fConversionCuts(NULL),
89  fClusterCutArray(NULL),
90  fCaloPhotonCuts(NULL),
91  fMesonCutArray(NULL),
92  fMesonCuts(NULL),
93  fHistoConvGammaPt(NULL),
94  fHistoMotherInvMassPt(NULL),
95  fHistoMotherMatchedInvMassPt(NULL),
96  fSparseMotherInvMassPtZM(NULL),
97  fHistoMotherBackInvMassPt(NULL),
98  fSparseMotherBackInvMassPtZM(NULL),
99  fHistoMotherMesonPtY(NULL),
100  fHistoMotherMesonPtAlpha(NULL),
101  fHistoMotherMesonPtOpenAngle(NULL),
102  fHistoMotherMesonConvPhotonEtaPhi(NULL),
103  fHistoClusGammaPt(NULL),
104  fHistoClusGammaE(NULL),
105  fHistoClusOverlapHeadersGammaPt(NULL),
106  fHistoClusAllHeadersGammaPt(NULL),
107  fHistoClusRejectedHeadersGammaPt(NULL),
108  fHistoMCHeaders(NULL),
109  fHistoMCMesonPt(NULL),
110  fHistoMCMesonWOWeightPt(NULL),
111  fHistoMCMesonWOEvtWeightPt(NULL),
112  fHistoMCMesonInAccPt(NULL),
113  fHistoMCMesonWOWeightInAccPt(NULL),
114  fHistoMCMesonWOEvtWeightInAccPt(NULL),
115  fHistoMCMesonPtY(NULL),
116  fHistoMCMesonPtAlpha(NULL),
117  fHistoMCMesonPtJetPt(NULL),
118  fHistoTrueMesonInvMassPt(NULL),
119  fHistoTrueMesonMatchedInvMassPt(NULL),
120  fHistoTrueMesonCaloPhotonInvMassPt(NULL),
121  fHistoTrueMesonCaloConvertedPhotonInvMassPt(NULL),
122  fHistoTrueMesonCaloMixedPhotonConvPhotonInvMassPt(NULL),
123  fHistoTrueMesonCaloConvertedPhotonMatchedInvMassPt(NULL),
124  fHistoTrueMesonCaloElectronInvMassPt(NULL),
125  fHistoTrueMesonCaloMergedClusterInvMassPt(NULL),
126  fHistoTrueMesonCaloMergedClusterPartConvInvMassPt(NULL),
127  fHistoTruePrimaryMesonInvMassPt(NULL),
128  fHistoTruePrimaryMesonW0WeightingInvMassPt(NULL),
129  fProfileTruePrimaryMesonWeightsInvMassPt(NULL),
130  fHistoTruePrimaryMesonMCPtResolPt(NULL),
131  fHistoTrueMotherMesonConvPhotonEtaPhi(NULL),
132  fHistoTrueBckGGInvMassPt(NULL),
133  fHistoTrueBckFullMesonContainedInOneClusterInvMassPt(NULL),
134  fHistoTrueBckAsymEClustersInvMassPt(NULL),
135  fHistoTrueBckContInvMassPt(NULL),
136  fHistoTrueMesonPtY(NULL),
137  fHistoTrueMesonPtAlpha(NULL),
138  fHistoTrueMesonPtOpenAngle(NULL),
139  fHistoTrueConvGammaPt(NULL),
140  fHistoTruePrimaryConvGammaPt(NULL),
141  fHistoTrueClusGammaPt(NULL),
142  fHistoTrueClusConvGammaPt(NULL),
143  fHistoTrueClusConvGammaFullyPt(NULL),
144  fHistoTrueNLabelsInClusPt(NULL),
145  fHistoTruePrimaryClusGammaPt(NULL),
146  fHistoTruePrimaryClusConvGammaPt(NULL),
147  fVectorRecTrueMesons(0),
148  fHistoDoubleCountTrueMesonInvMassPt(NULL),
149  fHistoDoubleCountTrueConvGammaRPt(NULL),
150  fHistoDoubleCountTrueClusterGammaPt(NULL),
151  fVectorDoubleCountTrueMesons(0),
152  fVectorDoubleCountTrueConvGammas(0),
153  fVectorDoubleCountTrueClusterGammas(0),
154  fHistoMultipleCountTrueMeson(NULL),
155  fHistoMultipleCountTrueConvGamma(NULL),
156  fHistoMultipleCountTrueClusterGamma(NULL),
157  fMapMultipleCountTrueMesons(),
158  fMapMultipleCountTrueConvGammas(),
159  fMapMultipleCountTrueClusterGammas(),
160  fHistoNEvents(NULL),
161  fHistoNEventsWOWeight(NULL),
162  fHistoNGoodESDTracks(NULL),
163  fHistoVertexZ(NULL),
164  fHistoVertexX(NULL),
165  fHistoVertexY(NULL),
166  fHistoNGammaConvCandidates(NULL),
167  fHistoNGammaCaloCandidates(NULL),
168  fHistoSPDClusterTrackletBackground(NULL),
169  fHistoNV0Tracks(NULL),
170  fProfileEtaShift(NULL),
171  fProfileJetJetXSection(NULL),
172  fHistoJetJetNTrials(NULL),
173  fEventPlaneAngle(-100),
174  fRandom(0),
175  fNGammaCandidates(0),
176  fUnsmearedPx(NULL),
177  fUnsmearedPy(NULL),
178  fUnsmearedPz(NULL),
179  fUnsmearedE(NULL),
180  fMCEventPos(NULL),
181  fMCEventNeg(NULL),
182  fESDArrayPos(NULL),
183  fESDArrayNeg(NULL),
184  fnCuts(0),
185  fiCut(0),
186  fMoveParticleAccordingToVertex(kTRUE),
187  fIsHeavyIon(0),
188  fDoLightOutput(kFALSE),
189  fMesonRecoMode(-1),
190  fMesonType(-1),
191  fMesonPDG(0),
192  fMesonInvMassMin(0),
193  fMesonInvMassMax(0),
194  fMesonInvMassNBins(0),
195  fMesonInvMassWindow(NULL),
196  fDoMesonQA(0),
197  fDoPhotonQA(0),
198  fDoClusterQA(0),
199  fIsFromDesiredHeader(kTRUE),
200  fIsOverlappingWithOtherHeader(kFALSE),
201  fIsMC(0),
202  fDoTHnSparse(kTRUE),
203  fSetPlotHistsExtQA(kFALSE),
204  fWeightJetJetMC(1),
205  fDoConvGammaShowerShapeTree(kFALSE),
206  fEnableSortForClusMC(kFALSE),
207  fDoPrimaryTrackMatching(kFALSE),
208  fDoInvMassShowerShapeTree(kFALSE),
209  tBrokenFiles(NULL),
210  fFileNameBroken(NULL),
211  fAllowOverlapHeaders(kTRUE)
212 {
213 
214 }
215 
216 //________________________________________________________________________
217 AliAnalysisTaskHeavyNeutralMesonToGG::AliAnalysisTaskHeavyNeutralMesonToGG(const char *name):
218  AliAnalysisTaskSE(name),
219  fV0Reader(NULL),
220  fV0ReaderName("V0ReaderV1"),
221  fCorrTaskSetting(""),
222  fBGHandler(NULL),
223  fBGHandlerRP(NULL),
224  fBGClusHandler(NULL),
225  fBGClusHandlerRP(NULL),
226  fInputEvent(NULL),
227  fMCEvent(NULL),
228  fCutFolder(NULL),
229  fESDList(NULL),
230  fBackList(NULL),
231  fMotherList(NULL),
232  fTrueList(NULL),
233  fMCList(NULL),
234  fOutputContainer(0),
235  fReaderGammas(NULL),
236  fGammaCandidates(NULL),
237  fClusterCandidates(NULL),
238  fEventCutArray(NULL),
239  fEventCuts(NULL),
240  fCutArray(NULL),
241  fConversionCuts(NULL),
242  fClusterCutArray(NULL),
243  fCaloPhotonCuts(NULL),
244  fMesonCutArray(NULL),
245  fMesonCuts(NULL),
246  fHistoConvGammaPt(NULL),
247  fHistoMotherInvMassPt(NULL),
248  fHistoMotherMatchedInvMassPt(NULL),
249  fSparseMotherInvMassPtZM(NULL),
250  fHistoMotherBackInvMassPt(NULL),
251  fSparseMotherBackInvMassPtZM(NULL),
252  fHistoMotherMesonPtY(NULL),
253  fHistoMotherMesonPtAlpha(NULL),
254  fHistoMotherMesonPtOpenAngle(NULL),
255  fHistoMotherMesonConvPhotonEtaPhi(NULL),
256  fHistoClusGammaPt(NULL),
257  fHistoClusGammaE(NULL),
258  fHistoClusOverlapHeadersGammaPt(NULL),
259  fHistoClusAllHeadersGammaPt(NULL),
260  fHistoClusRejectedHeadersGammaPt(NULL),
261  fHistoMCHeaders(NULL),
262  fHistoMCMesonPt(NULL),
263  fHistoMCMesonWOWeightPt(NULL),
264  fHistoMCMesonWOEvtWeightPt(NULL),
265  fHistoMCMesonInAccPt(NULL),
266  fHistoMCMesonWOWeightInAccPt(NULL),
267  fHistoMCMesonWOEvtWeightInAccPt(NULL),
268  fHistoMCMesonPtY(NULL),
269  fHistoMCMesonPtAlpha(NULL),
270  fHistoMCMesonPtJetPt(NULL),
271  fHistoTrueMesonInvMassPt(NULL),
272  fHistoTrueMesonMatchedInvMassPt(NULL),
273  fHistoTrueMesonCaloPhotonInvMassPt(NULL),
274  fHistoTrueMesonCaloConvertedPhotonInvMassPt(NULL),
275  fHistoTrueMesonCaloMixedPhotonConvPhotonInvMassPt(NULL),
276  fHistoTrueMesonCaloConvertedPhotonMatchedInvMassPt(NULL),
277  fHistoTrueMesonCaloElectronInvMassPt(NULL),
278  fHistoTrueMesonCaloMergedClusterInvMassPt(NULL),
279  fHistoTrueMesonCaloMergedClusterPartConvInvMassPt(NULL),
280  fHistoTruePrimaryMesonInvMassPt(NULL),
281  fHistoTruePrimaryMesonW0WeightingInvMassPt(NULL),
282  fProfileTruePrimaryMesonWeightsInvMassPt(NULL),
283  fHistoTruePrimaryMesonMCPtResolPt(NULL),
284  fHistoTrueMotherMesonConvPhotonEtaPhi(NULL),
285  fHistoTrueBckGGInvMassPt(NULL),
286  fHistoTrueBckFullMesonContainedInOneClusterInvMassPt(NULL),
287  fHistoTrueBckAsymEClustersInvMassPt(NULL),
288  fHistoTrueBckContInvMassPt(NULL),
289  fHistoTrueMesonPtY(NULL),
290  fHistoTrueMesonPtAlpha(NULL),
291  fHistoTrueMesonPtOpenAngle(NULL),
292  fHistoTrueConvGammaPt(NULL),
293  fHistoTruePrimaryConvGammaPt(NULL),
294  fHistoTrueClusGammaPt(NULL),
295  fHistoTrueClusConvGammaPt(NULL),
296  fHistoTrueClusConvGammaFullyPt(NULL),
297  fHistoTrueNLabelsInClusPt(NULL),
298  fHistoTruePrimaryClusGammaPt(NULL),
299  fHistoTruePrimaryClusConvGammaPt(NULL),
300  fVectorRecTrueMesons(0),
301  fHistoDoubleCountTrueMesonInvMassPt(NULL),
302  fHistoDoubleCountTrueConvGammaRPt(NULL),
303  fHistoDoubleCountTrueClusterGammaPt(NULL),
304  fVectorDoubleCountTrueMesons(0),
305  fVectorDoubleCountTrueConvGammas(0),
306  fVectorDoubleCountTrueClusterGammas(0),
307  fHistoMultipleCountTrueMeson(NULL),
308  fHistoMultipleCountTrueConvGamma(NULL),
309  fHistoMultipleCountTrueClusterGamma(NULL),
310  fMapMultipleCountTrueMesons(),
311  fMapMultipleCountTrueConvGammas(),
312  fMapMultipleCountTrueClusterGammas(),
313  fHistoNEvents(NULL),
314  fHistoNEventsWOWeight(NULL),
315  fHistoNGoodESDTracks(NULL),
316  fHistoVertexZ(NULL),
317  fHistoVertexX(NULL),
318  fHistoVertexY(NULL),
319  fHistoNGammaConvCandidates(NULL),
320  fHistoNGammaCaloCandidates(NULL),
321  fHistoSPDClusterTrackletBackground(NULL),
322  fHistoNV0Tracks(NULL),
323  fProfileEtaShift(NULL),
324  fProfileJetJetXSection(NULL),
325  fHistoJetJetNTrials(NULL),
326  fEventPlaneAngle(-100),
327  fRandom(0),
328  fNGammaCandidates(0),
329  fUnsmearedPx(NULL),
330  fUnsmearedPy(NULL),
331  fUnsmearedPz(NULL),
332  fUnsmearedE(NULL),
333  fMCEventPos(NULL),
334  fMCEventNeg(NULL),
335  fESDArrayPos(NULL),
336  fESDArrayNeg(NULL),
337  fnCuts(0),
338  fiCut(0),
339  fMoveParticleAccordingToVertex(kTRUE),
340  fIsHeavyIon(0),
341  fDoLightOutput(kFALSE),
342  fMesonRecoMode(-1),
343  fMesonType(-1),
344  fMesonPDG(0),
345  fMesonInvMassMin(0),
346  fMesonInvMassMax(0),
347  fMesonInvMassNBins(0),
348  fMesonInvMassWindow(NULL),
349  fDoMesonQA(0),
350  fDoPhotonQA(0),
351  fDoClusterQA(0),
352  fIsFromDesiredHeader(kTRUE),
353  fIsOverlappingWithOtherHeader(kFALSE),
354  fIsMC(0),
355  fDoTHnSparse(kTRUE),
356  fSetPlotHistsExtQA(kFALSE),
357  fWeightJetJetMC(1),
358  fDoConvGammaShowerShapeTree(kFALSE),
359  fEnableSortForClusMC(kFALSE),
360  fDoPrimaryTrackMatching(kFALSE),
361  fDoInvMassShowerShapeTree(kFALSE),
362  tBrokenFiles(NULL),
363  fFileNameBroken(NULL),
364  fAllowOverlapHeaders(kTRUE)
365 {
366  // Define output slots here
367  DefineOutput(1, TList::Class());
368 }
369 
370 AliAnalysisTaskHeavyNeutralMesonToGG::~AliAnalysisTaskHeavyNeutralMesonToGG()
371 {
372  if(fGammaCandidates){
373  delete fGammaCandidates;
374  fGammaCandidates = 0x0;
375  }
376  if(fClusterCandidates){
377  delete fClusterCandidates;
378  fClusterCandidates = 0x0;
379  }
380  if(fBGHandler){
381  delete[] fBGHandler;
382  fBGHandler = 0x0;
383  }
384  if(fBGHandlerRP){
385  delete[] fBGHandlerRP;
386  fBGHandlerRP = 0x0;
387  }
388  if(fBGClusHandler){
389  delete[] fBGClusHandler;
390  fBGClusHandler = 0x0;
391  }
392  if(fBGClusHandlerRP){
393  delete[] fBGClusHandlerRP;
394  fBGClusHandlerRP = 0x0;
395  }
396 }
397 //___________________________________________________________
398 void AliAnalysisTaskHeavyNeutralMesonToGG::InitBack(){
399 
400  const Int_t nDim = 4;
401  Int_t nBins[nDim] = {800,300,7,4};
402  Double_t xMin[nDim] = {0,0, 0,0};
403  Double_t xMax[nDim] = {0.8,30,7,4};
404 
405  if(fDoTHnSparse){
406  fSparseMotherInvMassPtZM = new THnSparseF*[fnCuts];
407  fSparseMotherBackInvMassPtZM = new THnSparseF*[fnCuts];
408  }
409 
410  fBGHandler = new AliGammaConversionAODBGHandler*[fnCuts];
411  fBGHandlerRP = new AliConversionAODBGHandlerRP*[fnCuts];
412 
413  fBGClusHandler = new AliGammaConversionAODBGHandler*[fnCuts];
414  fBGClusHandlerRP = new AliConversionAODBGHandlerRP*[fnCuts];
415 
416  for(Int_t iCut = 0; iCut<fnCuts;iCut++){
417  if (((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->DoBGCalculation()){
418  TString cutstringEvent = ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber();
419  TString cutstringConvGamma = "";
420  TString cutstringCaloGamma = "";
421  if (fMesonRecoMode < 2) cutstringConvGamma = ((AliConversionPhotonCuts*)fCutArray->At(iCut))->GetCutNumber();
422  if (fMesonRecoMode > 0) cutstringCaloGamma = ((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutNumber();
423  TString cutstringMeson = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
424 
425  TString fullCutString = "";
426  if (fMesonRecoMode == 0) fullCutString = Form("%i_%s_%s_%s",fMesonRecoMode, cutstringEvent.Data(), cutstringConvGamma.Data(), cutstringMeson.Data());
427  else if (fMesonRecoMode == 1) fullCutString = Form("%i_%s_%s_%s_%s",fMesonRecoMode, cutstringEvent.Data(),cutstringConvGamma.Data(),cutstringCaloGamma.Data(), cutstringMeson.Data());
428  else if (fMesonRecoMode == 2) fullCutString = Form("%i_%s_%s_%s",fMesonRecoMode, cutstringEvent.Data(), cutstringCaloGamma.Data(), cutstringMeson.Data());
429 
430  Int_t collisionSystem = atoi((TString)(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber())(0,1));
431  Int_t centMin = atoi((TString)(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber())(1,1));
432  Int_t centMax = atoi((TString)(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber())(2,1));
433 
434  if( collisionSystem == 1 || collisionSystem == 2 ||
435  collisionSystem == 5 || collisionSystem == 8 ||
436  collisionSystem == 9){
437  centMin = centMin*10;
438  centMax = centMax*10;
439  if(centMax ==0 && centMax!=centMin) centMax=100;
440  }else if(collisionSystem == 3 || collisionSystem == 6){
441  centMin = centMin*5;
442  centMax = centMax*5;
443  }else if(collisionSystem == 4 || collisionSystem == 7){
444  centMin = ((centMin*5)+45);
445  centMax = ((centMax*5)+45);
446  }
447 
448  if(fDoTHnSparse){
449  fBackList[iCut] = new TList();
450  fBackList[iCut]->SetName(Form("%s Back histograms",fullCutString.Data()));
451  fBackList[iCut]->SetOwner(kTRUE);
452  fCutFolder[iCut]->Add(fBackList[iCut]);
453 
454  fSparseMotherBackInvMassPtZM[iCut] = new THnSparseF("Back_Back_InvMass_Pt_z_m", "Back_Back_InvMass_Pt_z_m",nDim,nBins,xMin,xMax);
455  fBackList[iCut]->Add(fSparseMotherBackInvMassPtZM[iCut]);
456 
457  fMotherList[iCut] = new TList();
458  fMotherList[iCut]->SetName(Form("%s Mother histograms",fullCutString.Data()));
459  fMotherList[iCut]->SetOwner(kTRUE);
460  fCutFolder[iCut]->Add(fMotherList[iCut]);
461 
462  fSparseMotherInvMassPtZM[iCut] = new THnSparseF("Back_Mother_InvMass_Pt_z_m", "Back_Mother_InvMass_Pt_z_m",nDim,nBins,xMin,xMax);
463  fMotherList[iCut]->Add(fSparseMotherInvMassPtZM[iCut]);
464  }
465 
466  if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->BackgroundHandlerType() == 0){
467  fBGHandler[iCut] = new AliGammaConversionAODBGHandler( collisionSystem,centMin,centMax,
468  ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetNumberOfBGEvents(),
469  ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseTrackMultiplicity(),
470  2,8,5);
471  fBGClusHandler[iCut] = new AliGammaConversionAODBGHandler( collisionSystem,centMin,centMax,
472  ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetNumberOfBGEvents(),
473  ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseTrackMultiplicity(),
474  2,8,5);
475  fBGHandlerRP[iCut] = NULL;
476  }else{
477  fBGHandlerRP[iCut] = new AliConversionAODBGHandlerRP(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsHeavyIon(),
478  ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity(),
479  ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetNumberOfBGEvents());
480  fBGClusHandlerRP[iCut] = new AliConversionAODBGHandlerRP( ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsHeavyIon(),
481  ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity(),
482  ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetNumberOfBGEvents());
483  fBGHandler[iCut] = NULL;
484  }
485  }
486  }
487 }
488 
489 //________________________________________________________________________
490 void AliAnalysisTaskHeavyNeutralMesonToGG::UserCreateOutputObjects(){
491 
492  fV0Reader = (AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask(fV0ReaderName.Data());
493  if(!fV0Reader){printf("Error: No V0 Reader");return;}// GetV0Reader
494  if (fIsMC == 2){
495  fDoClusterQA = 0;
496  fDoTHnSparse = kFALSE;
497  } else if (fIsMC == 3){
498  fDoTHnSparse = kFALSE;
499  }
500  // Create histograms
501  if(fOutputContainer != NULL){
502  delete fOutputContainer;
503  fOutputContainer = NULL;
504  }
505  if(fOutputContainer == NULL){
506  fOutputContainer = new TList();
507  fOutputContainer->SetOwner(kTRUE);
508  }
509 
510  // Array of current cut's gammas
511  fGammaCandidates = new TList();
512  fClusterCandidates = new TList();
513  fClusterCandidates->SetOwner(kTRUE);
514 
515  fCutFolder = new TList*[fnCuts];
516  fESDList = new TList*[fnCuts];
517 
518  if(fDoTHnSparse){
519  fBackList = new TList*[fnCuts];
520  fMotherList = new TList*[fnCuts];
521  }
522 
523  fHistoNEvents = new TH1F*[fnCuts];
524  if (fIsMC > 1){
525  fHistoNEventsWOWeight = new TH1F*[fnCuts];
526  }
527  if (fIsMC == 2){
528  fProfileJetJetXSection = new TProfile*[fnCuts];
529  fHistoJetJetNTrials = new TH1F*[fnCuts];
530  }
531  fHistoNGoodESDTracks = new TH1F*[fnCuts];
532  fHistoVertexZ = new TH1F*[fnCuts];
533  if(!fDoLightOutput){
534  fHistoVertexX = new TH1F*[fnCuts];
535  fHistoVertexY = new TH1F*[fnCuts];
536  }
537  if (fMesonRecoMode < 2) fHistoNGammaConvCandidates = new TH1F*[fnCuts];
538  if(fIsHeavyIon==2) fProfileEtaShift = new TProfile*[fnCuts];
539  if(!fDoLightOutput){
540  fHistoSPDClusterTrackletBackground = new TH2F*[fnCuts];
541  fHistoNV0Tracks = new TH1F*[fnCuts];
542  if (fMesonRecoMode < 2)
543  fHistoConvGammaPt = new TH1F*[fnCuts];
544  }
545 
546  if (fMesonRecoMode > 0) {
547  fHistoNGammaCaloCandidates = new TH1F*[fnCuts];
548  if (!fDoLightOutput ){
549  fHistoClusGammaPt = new TH1F*[fnCuts];
550  fHistoClusGammaE = new TH1F*[fnCuts];
551  if (fIsMC > 0){
552  fHistoClusOverlapHeadersGammaPt = new TH1F*[fnCuts];
553  fHistoClusAllHeadersGammaPt = new TH1F*[fnCuts];
554  fHistoClusRejectedHeadersGammaPt = new TH1F*[fnCuts];
555  }
556  }
557  }
558 
559  fHistoMotherInvMassPt = new TH2F*[fnCuts];
560  fHistoMotherBackInvMassPt = new TH2F*[fnCuts];
561  if(!fDoLightOutput && fMesonRecoMode == 1){
562  fHistoMotherMatchedInvMassPt = new TH2F*[fnCuts];
563  }
564  if (fDoMesonQA > 0){
565  fHistoMotherMesonPtY = new TH2F*[fnCuts];
566  fHistoMotherMesonPtAlpha = new TH2F*[fnCuts];
567  fHistoMotherMesonPtOpenAngle = new TH2F*[fnCuts];
568  fHistoMotherMesonConvPhotonEtaPhi = new TH2F*[fnCuts];
569  }
570 
571  fMesonInvMassWindow = new Double_t[2];
572  if (fMesonType < 0 || fMesonType > 2){
573  if(!fV0Reader){printf("Error: No V0 Reader");return;}// GetV0Reader
574  } else if (fMesonType == 0){ // pi0 case 134.9770 ± 0.0005 MeV
575  fMesonPDG = 111;
576  fMesonInvMassMin = 0.;
577  fMesonInvMassMax = 0.400;
578  fMesonInvMassNBins = 400;
579  fMesonInvMassWindow[0] = 0.05;
580  fMesonInvMassWindow[1] = 0.17;
581  } else if (fMesonType == 1){ // eta case 547.862 ± 0.017 MeV
582  fMesonPDG = 221;
583  fMesonInvMassMin = 0.300;
584  fMesonInvMassMax = 0.800;
585  fMesonInvMassNBins = 500;
586  fMesonInvMassWindow[0] = 0.45;
587  fMesonInvMassWindow[1] = 0.65;
588  } else if (fMesonType == 2){ // eta' case 957.78 ± 0.06 MeV
589  fMesonPDG = 331;
590  fMesonInvMassMin = 0.700;
591  fMesonInvMassMax = 1.200;
592  fMesonInvMassNBins = 500;
593  fMesonInvMassWindow[0] = 0.85;
594  fMesonInvMassWindow[1] = 1.05;
595  }
596  // set common binning in pT for mesons and photons
597  Int_t nBinsPt = 200;
598  Float_t binWidthPt = 0.1;
599  Float_t minPt = 0;
600  Float_t maxPt = 20;
601  Int_t nBinsQAPt = 170;
602  Float_t maxQAPt = 20;
603  Int_t nBinsClusterPt = 500;
604  Float_t minClusterPt = 0;
605  Float_t maxClusterPt = 50;
606  Double_t *arrPtBinning = new Double_t[1200];
607  Double_t *arrQAPtBinning = new Double_t[1200];
608  Double_t *arrClusPtBinning = new Double_t[1200];
609  // Set special pt binning for pp 8TeV
610  if (((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEnergyEnum() == AliConvEventCuts::k8TeV ){
611  nBinsQAPt = 190;
612  maxQAPt = 40;
613  for(Int_t i=0; i<nBinsQAPt+1;i++){
614  if(i<60) arrQAPtBinning[i] = 0.05*i;
615  else if(i<130) arrQAPtBinning[i] = 3.+0.1*(i-60);
616  else if(i<170) arrQAPtBinning[i] = 10.+0.25*(i-130);
617  else if(i<190) arrQAPtBinning[i] = 20.+1.0*(i-170);
618  else arrQAPtBinning[i] = maxQAPt;
619  }
620  nBinsPt = 400;
621  minPt = 0;
622  maxPt = 40;
623  for(Int_t i=0; i<nBinsPt+1;i++){
624  arrPtBinning[i] = ((maxPt-minPt)/nBinsPt)*i;
625  }
626  nBinsClusterPt = 800;
627  minClusterPt = 0;
628  maxClusterPt = 80;
629  for(Int_t i=0; i<nBinsPt+1;i++){
630  arrClusPtBinning[i] = ((maxClusterPt-minClusterPt)/nBinsClusterPt)*i;
631  }
632  // Set special pt binning for pPb 5TeV
633  } else if ( ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEnergyEnum() == AliConvEventCuts::kpPb5TeV ||
634  ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEnergyEnum() == AliConvEventCuts::kpPb5TeVR2 ||
635  ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEnergyEnum() == AliConvEventCuts::k5TeV ){
636  binWidthPt = 0.05;
637  nBinsPt = 205;
638  minPt = 0;
639  maxPt = 60;
640  for(Int_t i=0; i<nBinsPt+1;i++){
641  if (i < 1) arrPtBinning[i] = 0.3*i;
642  else if(i<55) arrPtBinning[i] = 0.3+0.05*(i-1);
643  else if(i<125) arrPtBinning[i] = 3.+0.1*(i-55);
644  else if(i<165) arrPtBinning[i] = 10.+0.25*(i-125);
645  else if(i<205) arrPtBinning[i] = 20.+1.0*(i-165);
646  else arrPtBinning[i] = maxPt;
647  }
648  nBinsQAPt = 210;
649  maxQAPt = 60;
650  for(Int_t i=0; i<nBinsQAPt+1;i++){
651  if(i<60) arrQAPtBinning[i] = 0.05*i;
652  else if(i<130) arrQAPtBinning[i] = 3.+0.1*(i-60);
653  else if(i<170) arrQAPtBinning[i] = 10.+0.25*(i-130);
654  else if(i<210) arrQAPtBinning[i] = 20.+1.0*(i-170);
655  else arrQAPtBinning[i] = maxQAPt;
656  }
657  nBinsClusterPt = 301;
658  minClusterPt = 0;
659  maxClusterPt = 100;
660  for(Int_t i=0; i<nBinsClusterPt+1;i++){
661  if (i < 1) arrClusPtBinning[i] = 0.3*i;
662  else if(i<55) arrClusPtBinning[i] = 0.3+0.05*(i-1);
663  else if(i<125) arrClusPtBinning[i] = 3.+0.1*(i-55);
664  else if(i<155) arrClusPtBinning[i] = 10.+0.2*(i-125);
665  else if(i<211) arrClusPtBinning[i] = 16.+0.25*(i-155);
666  else if(i<251) arrClusPtBinning[i] = 30.+0.5*(i-211);
667  else if(i<301) arrClusPtBinning[i] = 50.+1.0*(i-251);
668  else arrClusPtBinning[i] = maxClusterPt;
669  }
670  // Set special pt binning for pp 13TeV, pPb 8TeV
671  } else if ( ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEnergyEnum() == AliConvEventCuts::k13TeV ||
672  ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEnergyEnum() == AliConvEventCuts::k13TeVLowB ||
673  ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEnergyEnum() == AliConvEventCuts::kpPb8TeV ){
674  nBinsPt = 285;
675  minPt = 0;
676  maxPt = 100;
677  binWidthPt = 0.05;
678  for(Int_t i=0; i<nBinsPt+1;i++){
679  if (i < 1) arrPtBinning[i] = 0.3*i;
680  else if(i<55) arrPtBinning[i] = 0.3+0.05*(i-1);
681  else if(i<125) arrPtBinning[i] = 3.+0.1*(i-55);
682  else if(i<185) arrPtBinning[i] = 10.+0.25*(i-125);
683  else if(i<235) arrPtBinning[i] = 25.+0.5*(i-185);
684  else if(i<285) arrPtBinning[i] = 50.+1.0*(i-235);
685  else arrPtBinning[i] = maxPt;
686  }
687  nBinsQAPt = 270;
688  maxQAPt = 100;
689  for(Int_t i=0; i<nBinsQAPt+1;i++){
690  if(i<60) arrQAPtBinning[i] = 0.05*i;
691  else if(i<130) arrQAPtBinning[i] = 3.+0.1*(i-60);
692  else if(i<170) arrQAPtBinning[i] = 10.+0.25*(i-130);
693  else if(i<210) arrQAPtBinning[i] = 20.+0.5*(i-170);
694  else if(i<270) arrQAPtBinning[i] = 40.+1.0*(i-210);
695  else arrQAPtBinning[i] = maxQAPt;
696  }
697  nBinsClusterPt = 301;
698  minClusterPt = 0;
699  maxClusterPt = 100;
700  for(Int_t i=0; i<nBinsClusterPt+1;i++){
701  if (i < 1) arrClusPtBinning[i] = 0.3*i;
702  else if(i<55) arrClusPtBinning[i] = 0.3+0.05*(i-1);
703  else if(i<125) arrClusPtBinning[i] = 3.+0.1*(i-55);
704  else if(i<155) arrClusPtBinning[i] = 10.+0.2*(i-125);
705  else if(i<211) arrClusPtBinning[i] = 16.+0.25*(i-155);
706  else if(i<251) arrClusPtBinning[i] = 30.+0.5*(i-211);
707  else if(i<301) arrClusPtBinning[i] = 50.+1.0*(i-251);
708  else arrClusPtBinning[i] = maxClusterPt;
709  }
710  } else if (((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEnergyEnum() == AliConvEventCuts::kXeXe5440GeV ){
711  nBinsPt = 90;
712  minPt = 0;
713  maxPt = 20;
714  for(Int_t i=0; i<nBinsPt+1;i++){
715  if (i < 1) arrPtBinning[i] = 0.3*i;
716  else if(i<58) arrPtBinning[i] = 0.3+0.1*(i-1);
717  else if(i<82) arrPtBinning[i] = 6.+0.25*(i-58);
718  else if(i<90) arrPtBinning[i] = 12.+1.0*(i-82);
719  else arrPtBinning[i] = maxPt;
720  }
721  nBinsQAPt = 92;
722  maxQAPt = 20;
723  for(Int_t i=0; i<nBinsQAPt+1;i++){
724  if(i<60) arrQAPtBinning[i] = 0.1*i;
725  else if(i<84) arrQAPtBinning[i] = 6.+0.25*(i-60);
726  else if(i<92) arrQAPtBinning[i] = 12.+1.0*(i-84);
727  else arrQAPtBinning[i] = maxQAPt;
728  }
729  nBinsClusterPt = 148;
730  minClusterPt = 0;
731  maxClusterPt = 40;
732  for(Int_t i=0; i<nBinsClusterPt+1;i++){
733  if (i < 1) arrClusPtBinning[i] = 0.3*i;
734  else if(i<98) arrClusPtBinning[i] = 0.3+0.1*(i-1);
735  else if(i<123) arrClusPtBinning[i] = 10.+0.2*(i-98);
736  else if(i<148) arrClusPtBinning[i] = 15.+1.0*(i-123);
737  else arrClusPtBinning[i] = maxClusterPt;
738  }
739  } else if (((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEnergyEnum() == AliConvEventCuts::kPbPb5TeV ){
740  nBinsPt = 90;
741  minPt = 0;
742  maxPt = 20;
743  for(Int_t i=0; i<nBinsPt+1;i++){
744  if (i < 1) arrPtBinning[i] = 0.3*i;
745  else if(i<58) arrPtBinning[i] = 0.3+0.1*(i-1);
746  else if(i<82) arrPtBinning[i] = 6.+0.25*(i-58);
747  else if(i<90) arrPtBinning[i] = 12.+1.0*(i-82);
748  else arrPtBinning[i] = maxPt;
749  }
750  nBinsQAPt = 92;
751  maxQAPt = 20;
752  for(Int_t i=0; i<nBinsQAPt+1;i++){
753  if(i<60) arrQAPtBinning[i] = 0.1*i;
754  else if(i<84) arrQAPtBinning[i] = 6.+0.25*(i-60);
755  else if(i<92) arrQAPtBinning[i] = 12.+1.0*(i-84);
756  else arrQAPtBinning[i] = maxQAPt;
757  }
758  nBinsClusterPt = 148;
759  minClusterPt = 0;
760  maxClusterPt = 40;
761  for(Int_t i=0; i<nBinsClusterPt+1;i++){
762  if (i < 1) arrClusPtBinning[i] = 0.3*i;
763  else if(i<98) arrClusPtBinning[i] = 0.3+0.1*(i-1);
764  else if(i<123) arrClusPtBinning[i] = 10.+0.2*(i-98);
765  else if(i<148) arrClusPtBinning[i] = 15.+1.0*(i-123);
766  else arrClusPtBinning[i] = maxClusterPt;
767  }
768  // default binning
769  } else {
770  for(Int_t i=0; i<nBinsPt+1;i++){
771  arrPtBinning[i] = ((maxPt-minPt)/nBinsPt)*i;
772  }
773  for(Int_t i=0; i<nBinsClusterPt+1;i++){
774  arrClusPtBinning[i] = ((maxClusterPt-minClusterPt)/nBinsClusterPt)*i;
775  }
776  for(Int_t i=0; i<nBinsQAPt+1;i++){
777  if(i<60) arrQAPtBinning[i] = 0.05*i;
778  else if(i<130) arrQAPtBinning[i] = 3.+0.1*(i-60);
779  else if(i<170) arrQAPtBinning[i] = 10.+0.25*(i-130);
780  else arrQAPtBinning[i] = maxQAPt;
781  }
782  }
783 
784  for(Int_t iCut = 0; iCut<fnCuts;iCut++){
785  TString cutstringEvent = ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber();
786  TString cutstringConvGamma = "";
787  TString cutstringCaloGamma = "";
788  if (fMesonRecoMode < 2) cutstringConvGamma = ((AliConversionPhotonCuts*)fCutArray->At(iCut))->GetCutNumber();
789  if (fMesonRecoMode > 0) cutstringCaloGamma = ((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutNumber();
790  TString cutstringMeson = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
791 
792  TString fullCutString = "";
793  if (fMesonRecoMode == 0) fullCutString = Form("%i_%s_%s_%s",fMesonRecoMode, cutstringEvent.Data(), cutstringConvGamma.Data(), cutstringMeson.Data());
794  else if (fMesonRecoMode == 1) fullCutString = Form("%i_%s_%s_%s_%s",fMesonRecoMode, cutstringEvent.Data(),cutstringConvGamma.Data(),cutstringCaloGamma.Data(), cutstringMeson.Data());
795  else if (fMesonRecoMode == 2) fullCutString = Form("%i_%s_%s_%s",fMesonRecoMode, cutstringEvent.Data(), cutstringCaloGamma.Data(), cutstringMeson.Data());
796 
797  fCutFolder[iCut] = new TList();
798  fCutFolder[iCut]->SetName(Form("Cut Number %s",fullCutString.Data()));
799  fCutFolder[iCut]->SetOwner(kTRUE);
800  fOutputContainer->Add(fCutFolder[iCut]);
801  fESDList[iCut] = new TList();
802  fESDList[iCut]->SetName(Form("%s ESD histograms",fullCutString.Data()));
803  fESDList[iCut]->SetOwner(kTRUE);
804  fCutFolder[iCut]->Add(fESDList[iCut]);
805 
806  fHistoNEvents[iCut] = new TH1F("NEvents", "NEvents", 14, -0.5, 13.5);
807  fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(1,"Accepted");
808  fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(2,"Centrality");
809  fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(3,"Miss. MC or inc. ev.");
810  if (((AliConvEventCuts*)fEventCutArray->At(iCut))->IsSpecialTrigger() > 1 ){
811  TString TriggerNames = "Not Trigger: ";
812  TriggerNames = TriggerNames+ ( (AliConvEventCuts*)fEventCutArray->At(iCut))->GetSpecialTriggerName();
813  fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(4,TriggerNames.Data());
814  }else {
815  fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(4,"Trigger");
816  }
817  fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(5,"Vertex Z");
818  fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(6,"Cont. Vertex");
819  fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(7,"Pile-Up");
820  fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(8,"no SDD");
821  fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(9,"no V0AND");
822  fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(10,"EMCAL problem");
823  fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(11,"rejectedForJetJetMC");
824  fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(12,"SPD hits vs tracklet");
825  fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(13,"Out-of-Bunch pileup Past-Future");
826  fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(14,"Pileup V0M-TPCout Tracks");
827  fESDList[iCut]->Add(fHistoNEvents[iCut]);
828 
829  if (fIsMC > 1){
830  fHistoNEventsWOWeight[iCut] = new TH1F("NEventsWOWeight", "NEventsWOWeight", 14, -0.5, 13.5);
831  fHistoNEventsWOWeight[iCut]->GetXaxis()->SetBinLabel(1,"Accepted");
832  fHistoNEventsWOWeight[iCut]->GetXaxis()->SetBinLabel(2,"Centrality");
833  fHistoNEventsWOWeight[iCut]->GetXaxis()->SetBinLabel(3,"Miss. MC or inc. ev.");
834  if (((AliConvEventCuts*)fEventCutArray->At(iCut))->IsSpecialTrigger() > 1 ){
835  TString TriggerNames = "Not Trigger: ";
836  TriggerNames = TriggerNames+ ( (AliConvEventCuts*)fEventCutArray->At(iCut))->GetSpecialTriggerName();
837  fHistoNEventsWOWeight[iCut]->GetXaxis()->SetBinLabel(4,TriggerNames.Data());
838  }else {
839  fHistoNEventsWOWeight[iCut]->GetXaxis()->SetBinLabel(4,"Trigger");
840  }
841  fHistoNEventsWOWeight[iCut]->GetXaxis()->SetBinLabel(5,"Vertex Z");
842  fHistoNEventsWOWeight[iCut]->GetXaxis()->SetBinLabel(6,"Cont. Vertex");
843  fHistoNEventsWOWeight[iCut]->GetXaxis()->SetBinLabel(7,"Pile-Up");
844  fHistoNEventsWOWeight[iCut]->GetXaxis()->SetBinLabel(8,"no SDD");
845  fHistoNEventsWOWeight[iCut]->GetXaxis()->SetBinLabel(9,"no V0AND");
846  fHistoNEventsWOWeight[iCut]->GetXaxis()->SetBinLabel(10,"EMCAL problem");
847  fHistoNEventsWOWeight[iCut]->GetXaxis()->SetBinLabel(11,"rejectedForJetJetMC");
848  fHistoNEventsWOWeight[iCut]->GetXaxis()->SetBinLabel(12,"SPD hits vs tracklet");
849  fHistoNEventsWOWeight[iCut]->GetXaxis()->SetBinLabel(13,"Out-of-Bunch pileup Past-Future");
850  fHistoNEventsWOWeight[iCut]->GetXaxis()->SetBinLabel(14,"Pileup V0M-TPCout Tracks");
851  fESDList[iCut]->Add(fHistoNEventsWOWeight[iCut]);
852  }
853 
854  if (fIsMC == 2){
855  fProfileJetJetXSection[iCut] = new TProfile("XSection", "XSection", 1, -0.5, 0.5);
856  fESDList[iCut]->Add(fProfileJetJetXSection[iCut]);
857  fHistoJetJetNTrials[iCut] = new TH1F("NTrials", "#sum{NTrials}", 1, 0, 1);
858  fHistoJetJetNTrials[iCut]->GetXaxis()->SetBinLabel(1,"#sum{NTrials}");
859  fESDList[iCut]->Add(fHistoJetJetNTrials[iCut]);
860  }
861 
862  if(fIsHeavyIon == 1)
863  fHistoNGoodESDTracks[iCut] = new TH1F("GoodESDTracks", "GoodESDTracks; # TPC tracks", 4000, 0, 4000);
864  else if(fIsHeavyIon == 2)
865  fHistoNGoodESDTracks[iCut] = new TH1F("GoodESDTracks", "GoodESDTracks; # TPC tracks", 400, 0, 400);
866  else
867  fHistoNGoodESDTracks[iCut] = new TH1F("GoodESDTracks", "GoodESDTracks; # TPC tracks", 200, 0, 200);
868  fESDList[iCut]->Add(fHistoNGoodESDTracks[iCut]);
869 
870  fHistoVertexZ[iCut] = new TH1F("VertexZ", "VertexZ", 200, -10, 10);
871  fESDList[iCut]->Add(fHistoVertexZ[iCut]);
872  if(!fDoLightOutput){
873  fHistoVertexX[iCut] = new TH1F("VertexX", "VertexX", 100, -5, 5);
874  fESDList[iCut]->Add(fHistoVertexX[iCut]);
875  fHistoVertexY[iCut] = new TH1F("VertexY", "VertexY", 100, -5, 5);
876  fESDList[iCut]->Add(fHistoVertexY[iCut]);
877  }
878 
879  if (fMesonRecoMode < 2){
880  if(fIsHeavyIon == 1)
881  fHistoNGammaConvCandidates[iCut] = new TH1F("GammaConvCandidates", "GammaConvCandidates; # accepted #gamma_{conv}", 100, 0, 100);
882  else if(fIsHeavyIon == 2)
883  fHistoNGammaConvCandidates[iCut] = new TH1F("GammaConvCandidates", "GammaConvCandidates; # accepted #gamma_{conv}", 50, 0, 50);
884  else
885  fHistoNGammaConvCandidates[iCut] = new TH1F("GammaConvCandidates", "GammaConvCandidates; # accepted #gamma_{conv}", 50, 0, 50);
886  fESDList[iCut]->Add(fHistoNGammaConvCandidates[iCut]);
887  }
888  if (fMesonRecoMode > 0){
889  if(fIsHeavyIon == 1)
890  fHistoNGammaCaloCandidates[iCut] = new TH1F("GammaCaloCandidates", "GammaCaloCandidates; # accepted #gamma_{conv}", 100, 0, 100);
891  else if(fIsHeavyIon == 2)
892  fHistoNGammaCaloCandidates[iCut] = new TH1F("GammaCaloCandidates", "GammaCaloCandidates; # accepted #gamma_{conv}", 50, 0, 50);
893  else
894  fHistoNGammaCaloCandidates[iCut] = new TH1F("GammaCaloCandidates", "GammaCaloCandidates; # accepted #gamma_{conv}", 50, 0, 50);
895  fESDList[iCut]->Add(fHistoNGammaCaloCandidates[iCut]);
896  }
897 
898  if(!fDoLightOutput){
899  fHistoSPDClusterTrackletBackground[iCut] = new TH2F("SPD tracklets vs SPD clusters", "SPD tracklets vs SPD clusters", 100, 0, 200, 250, 0, 1000);
900  fESDList[iCut]->Add(fHistoSPDClusterTrackletBackground[iCut]);
901 
902  if(fIsHeavyIon == 1)
903  fHistoNV0Tracks[iCut] = new TH1F("V0 Multiplicity", "V0 Multiplicity; VZERO amp [arb. units]", 30000, 0, 30000);
904  else if(fIsHeavyIon == 2)
905  fHistoNV0Tracks[iCut] = new TH1F("V0 Multiplicity", "V0 Multiplicity; VZERO amp [arb. units]", 2500, 0, 2500);
906  else
907  fHistoNV0Tracks[iCut] = new TH1F("V0 Multiplicity", "V0 Multiplicity; VZERO amp [arb. units]", 1500, 0, 1500);
908  fESDList[iCut]->Add(fHistoNV0Tracks[iCut]);
909 
910  if (fMesonRecoMode < 2){
911  fHistoConvGammaPt[iCut] = new TH1F("ESD_ConvGamma_Pt", "ESD_ConvGamma_Pt; p_{T,conv}(GeV/c)", (Int_t)((maxPt-minPt)/binWidthPt), minPt, maxPt);
912  fESDList[iCut]->Add(fHistoConvGammaPt[iCut]);
913  }
914  if (fMesonRecoMode > 0){
915  fHistoClusGammaPt[iCut] = new TH1F("ClusGamma_Pt", "ClusGamma_Pt; p_{T,clus} (GeV/c)", nBinsClusterPt, arrClusPtBinning);
916  fESDList[iCut]->Add(fHistoClusGammaPt[iCut]);
917  fHistoClusGammaE[iCut] = new TH1F("ClusGamma_E", "ClusGamma_E; E_{clus} (GeV)", nBinsClusterPt, arrClusPtBinning);
918  fESDList[iCut]->Add(fHistoClusGammaE[iCut]);
919  if (fIsMC > 0){
920  fHistoClusOverlapHeadersGammaPt[iCut] = new TH1F("ClusGammaOverlapHeaders_Pt", "ClusGammaOverlapHeaders_Pt; p_{T,clus} (GeV/c), selected header w/ overlap",
921  nBinsClusterPt, arrClusPtBinning);
922  fESDList[iCut]->Add(fHistoClusOverlapHeadersGammaPt[iCut]);
923  fHistoClusAllHeadersGammaPt[iCut] = new TH1F("ClusGammaAllHeaders_Pt", "ClusGammaAllHeaders_Pt; p_{T,clus} (GeV/c), all headers",
924  nBinsClusterPt, arrClusPtBinning);
925  fESDList[iCut]->Add(fHistoClusAllHeadersGammaPt[iCut]);
926  fHistoClusRejectedHeadersGammaPt[iCut] = new TH1F("ClusGammaRejectedHeaders_Pt", "ClusGammaRejectedHeaders_Pt; p_{T,clus} (GeV/c), rejected headers",
927  nBinsClusterPt, arrClusPtBinning);
928  fESDList[iCut]->Add(fHistoClusRejectedHeadersGammaPt[iCut]);
929  }
930  }
931  }
932 
933  if(fIsHeavyIon == 2){
934  fProfileEtaShift[iCut] = new TProfile("Eta Shift", "Eta Shift", 1, -0.5, 0.5);
935  fESDList[iCut]->Add(fProfileEtaShift[iCut]);
936  }
937 
938  if (fIsMC > 1){
939  fHistoNEvents[iCut]->Sumw2();
940  fHistoNGoodESDTracks[iCut]->Sumw2();
941  fHistoVertexZ[iCut]->Sumw2();
942  if (fHistoNGammaConvCandidates[iCut]) fHistoNGammaConvCandidates[iCut]->Sumw2();
943  if (fHistoNGammaCaloCandidates[iCut]) fHistoNGammaCaloCandidates[iCut]->Sumw2();
944  if(!fDoLightOutput){
945  fHistoVertexX[iCut]->Sumw2();
946  fHistoVertexY[iCut]->Sumw2();
947  fHistoSPDClusterTrackletBackground[iCut]->Sumw2();
948  fHistoNV0Tracks[iCut]->Sumw2();
949  if (fHistoConvGammaPt[iCut]) fHistoConvGammaPt[iCut]->Sumw2();
950  if (fHistoClusGammaPt[iCut]) fHistoClusGammaPt[iCut]->Sumw2();
951  if (fHistoClusGammaE[iCut]) fHistoClusGammaE[iCut]->Sumw2();
952  if (fHistoClusOverlapHeadersGammaPt[iCut]) fHistoClusOverlapHeadersGammaPt[iCut]->Sumw2();
953  if (fHistoClusAllHeadersGammaPt[iCut]) fHistoClusAllHeadersGammaPt[iCut]->Sumw2();
954  if (fHistoClusRejectedHeadersGammaPt[iCut]) fHistoClusRejectedHeadersGammaPt[iCut]->Sumw2();
955  }
956  }
957 
958 
959  fHistoMotherInvMassPt[iCut] = new TH2F("ESD_Mother_InvMass_Pt", "ESD_Mother_InvMass_Pt; M_{inv} (GeV/c^{2}); p_{T,pair} (GeV/c)",
960  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
961  fESDList[iCut]->Add(fHistoMotherInvMassPt[iCut]);
962 
963  if(!fDoLightOutput){
964  if (fMesonRecoMode == 1){
965  fHistoMotherMatchedInvMassPt[iCut] = new TH2F("ESD_MotherMatched_InvMass_Pt", "ESD_MotherMatched_InvMass_Pt; M_{inv} (GeV/c^{2}) matched conv e^{+/-}to cluster; p_{T,pair} (GeV/c)",
966  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
967  fESDList[iCut]->Add(fHistoMotherMatchedInvMassPt[iCut]);
968  }
969  }
970 
971  fHistoMotherBackInvMassPt[iCut] = new TH2F("ESD_Background_InvMass_Pt", "ESD_Background_InvMass_Pt; M_{inv, mxed}(GeV/c^{2}); p_{T,BG pair} (GeV/c)",
972  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
973  fESDList[iCut]->Add(fHistoMotherBackInvMassPt[iCut]);
974 
975  if (fIsMC > 1){
976  fHistoMotherInvMassPt[iCut]->Sumw2();
977  fHistoMotherBackInvMassPt[iCut]->Sumw2();
978  if (fHistoMotherMatchedInvMassPt[iCut]) fHistoMotherMatchedInvMassPt[iCut]->Sumw2();
979  }
980 
981  if (fDoMesonQA > 0 ){
982  fHistoMotherMesonPtY[iCut] = new TH2F("ESD_MotherMeson_Pt_Y", "ESD_MotherMeson_Pt_Y; p_{T, meson cand} (GeV/c); y_{meson cand}",
983  nBinsQAPt, arrQAPtBinning, 150, -1.5, 1.5);
984  fESDList[iCut]->Add(fHistoMotherMesonPtY[iCut]);
985  fHistoMotherMesonPtAlpha[iCut] = new TH2F("ESD_MotherMeson_Pt_Alpha", "ESD_MotherMeson_Pt_Alpha; p_{T, meson cand} (GeV/c); #alpha_{meson cand}",
986  nBinsQAPt, arrQAPtBinning, 200, -1, 1);
987  fESDList[iCut]->Add(fHistoMotherMesonPtAlpha[iCut]);
988  fHistoMotherMesonPtOpenAngle[iCut] = new TH2F("ESD_MotherMeson_Pt_OpenAngle", "ESD_MotherMeson_Pt_OpenAngle; p_{T, meson cand} (GeV/c); #theta_{meson cand}",
989  nBinsQAPt, arrQAPtBinning, 100, 0, 1);
990  fESDList[iCut]->Add(fHistoMotherMesonPtOpenAngle[iCut]);
991  fHistoMotherMesonConvPhotonEtaPhi[iCut] = new TH2F("ESD_MotherMesonConvPhoton_Eta_Phi", "ConvPhoton under meson peak; #phi_{#gamma_{conv}}(rad); #eta_{#gamma_{conv}}",
992  600, 0, 2*TMath::Pi(), 200, -1, 1);
993  fESDList[iCut]->Add(fHistoMotherMesonConvPhotonEtaPhi[iCut]);
994  if (fIsMC > 1){
995  fHistoMotherMesonPtY[iCut]->Sumw2();
996  fHistoMotherMesonPtAlpha[iCut]->Sumw2();
997  fHistoMotherMesonPtOpenAngle[iCut]->Sumw2();
998  fHistoMotherMesonConvPhotonEtaPhi[iCut]->Sumw2();
999  }
1000  }
1001  }
1002 
1003  InitBack(); // Init Background Handler
1004 
1005  if(fIsMC>0){
1006  // MC Histogramms
1007  fMCList = new TList*[fnCuts];
1008  // True Histogramms
1009  fTrueList = new TList*[fnCuts];
1010 
1011  if(!fDoLightOutput){
1012  fHistoMCHeaders = new TH1I*[fnCuts];
1013  if (fMesonRecoMode < 2){
1014  fHistoTrueConvGammaPt = new TH1F*[fnCuts];
1015  fHistoDoubleCountTrueConvGammaRPt = new TH2F*[fnCuts];
1016  fHistoMultipleCountTrueConvGamma = new TH1F*[fnCuts];
1017  fHistoTruePrimaryConvGammaPt = new TH1F*[fnCuts];
1018  }
1019 
1020  if (fMesonRecoMode > 0){
1021  fHistoTrueClusGammaPt = new TH1F*[fnCuts];
1022  fHistoTrueClusConvGammaPt = new TH1F*[fnCuts];
1023  fHistoTruePrimaryClusGammaPt = new TH1F*[fnCuts];
1024  fHistoTruePrimaryClusConvGammaPt = new TH1F*[fnCuts];
1025  fHistoDoubleCountTrueClusterGammaPt = new TH2F*[fnCuts];
1026  fHistoMultipleCountTrueClusterGamma = new TH1F*[fnCuts];
1027  fHistoTrueNLabelsInClusPt = new TH2F*[fnCuts];
1028  fHistoTrueClusConvGammaFullyPt = new TH1F*[fnCuts];
1029  }
1030  }
1031 
1032  fHistoMCMesonPt = new TH1F*[fnCuts];
1033  fHistoMCMesonWOWeightPt = new TH1F*[fnCuts];
1034  fHistoMCMesonInAccPt = new TH1F*[fnCuts];
1035  fHistoMCMesonWOWeightInAccPt = new TH1F*[fnCuts];
1036  if (fIsMC > 1){
1037  fHistoMCMesonWOEvtWeightPt = new TH1F*[fnCuts];
1038  fHistoMCMesonWOEvtWeightInAccPt = new TH1F*[fnCuts];
1039  }
1040 
1041  fHistoTrueMesonInvMassPt = new TH2F*[fnCuts];
1042  fHistoDoubleCountTrueMesonInvMassPt = new TH2F*[fnCuts];
1043  fHistoMultipleCountTrueMeson = new TH1F*[fnCuts];
1044  fHistoTruePrimaryMesonInvMassPt = new TH2F*[fnCuts];
1045  fHistoTruePrimaryMesonW0WeightingInvMassPt = new TH2F*[fnCuts];
1046  fProfileTruePrimaryMesonWeightsInvMassPt = new TProfile2D*[fnCuts];
1047  if (fMesonRecoMode == 1){
1048  fHistoTrueMesonMatchedInvMassPt = new TH2F*[fnCuts];
1049  }
1050 
1051  if (fDoMesonQA > 0){
1052  fHistoMCMesonPtY = new TH2F*[fnCuts];
1053  fHistoMCMesonPtAlpha = new TH2F*[fnCuts];
1054  if (fIsMC == 2){
1055  fHistoMCMesonPtJetPt = new TH2F*[fnCuts];
1056  }
1057 
1058  if (fIsMC < 2){
1059  if (fMesonRecoMode > 0){
1060  fHistoTrueMesonCaloPhotonInvMassPt = new TH2F*[fnCuts];
1061  fHistoTrueMesonCaloConvertedPhotonInvMassPt = new TH2F*[fnCuts];
1062  fHistoTrueMesonCaloElectronInvMassPt = new TH2F*[fnCuts];
1063  fHistoTrueMesonCaloMergedClusterInvMassPt = new TH2F*[fnCuts];
1064  fHistoTrueMesonCaloMergedClusterPartConvInvMassPt = new TH2F*[fnCuts];
1065  }
1066  if (fMesonRecoMode == 1){
1067  fHistoTrueMesonCaloConvertedPhotonMatchedInvMassPt = new TH2F*[fnCuts];
1068  fHistoTrueMotherMesonConvPhotonEtaPhi = new TH2F*[fnCuts];
1069  }
1070  if (fMesonRecoMode == 2){
1071  fHistoTrueMesonCaloMixedPhotonConvPhotonInvMassPt = new TH2F*[fnCuts];
1072  }
1073  fHistoTruePrimaryMesonMCPtResolPt = new TH2F*[fnCuts];
1074  }
1075  if(fDoMesonQA > 1){
1076  fHistoTrueBckGGInvMassPt = new TH2F*[fnCuts];
1077  fHistoTrueBckFullMesonContainedInOneClusterInvMassPt = new TH2F*[fnCuts];
1078  fHistoTrueBckAsymEClustersInvMassPt = new TH2F*[fnCuts];
1079  fHistoTrueBckContInvMassPt = new TH2F*[fnCuts];
1080  }
1081  fHistoTrueMesonPtY = new TH2F*[fnCuts];
1082  fHistoTrueMesonPtAlpha = new TH2F*[fnCuts];
1083  fHistoTrueMesonPtOpenAngle = new TH2F*[fnCuts];
1084  }
1085 
1086  for(Int_t iCut = 0; iCut<fnCuts;iCut++){
1087  TString cutstringEvent = ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber();
1088  TString cutstringConvGamma = "";
1089  TString cutstringCaloGamma = "";
1090  if (fMesonRecoMode < 2) cutstringConvGamma = ((AliConversionPhotonCuts*)fCutArray->At(iCut))->GetCutNumber();
1091  if (fMesonRecoMode > 0) cutstringCaloGamma = ((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutNumber();
1092  TString cutstringMeson = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
1093 
1094  TString fullCutString = "";
1095  if (fMesonRecoMode == 0) fullCutString = Form("%i_%s_%s_%s",fMesonRecoMode, cutstringEvent.Data(), cutstringConvGamma.Data(), cutstringMeson.Data());
1096  else if (fMesonRecoMode == 1) fullCutString = Form("%i_%s_%s_%s_%s",fMesonRecoMode, cutstringEvent.Data(),cutstringConvGamma.Data(),cutstringCaloGamma.Data(), cutstringMeson.Data());
1097  else if (fMesonRecoMode == 2) fullCutString = Form("%i_%s_%s_%s",fMesonRecoMode, cutstringEvent.Data(), cutstringCaloGamma.Data(), cutstringMeson.Data());
1098 
1099  fMCList[iCut] = new TList();
1100  fMCList[iCut]->SetName(Form("%s MC histograms",fullCutString.Data()));
1101  fMCList[iCut]->SetOwner(kTRUE);
1102  fCutFolder[iCut]->Add(fMCList[iCut]);
1103  if(!fDoLightOutput){
1104  fHistoMCHeaders[iCut] = new TH1I("MC_Headers", "MC_Headers", 20, 0, 20);
1105  fMCList[iCut]->Add(fHistoMCHeaders[iCut]);
1106  }
1107 
1108  fHistoMCMesonPt[iCut] = new TH1F("MC_Meson_Pt", "MC_Meson_Pt; p_{T} (GeV/c)",
1109  (Int_t)((maxPt-minPt)/binWidthPt), minPt, maxPt);
1110  fHistoMCMesonPt[iCut]->Sumw2();
1111  fMCList[iCut]->Add(fHistoMCMesonPt[iCut]);
1112  fHistoMCMesonWOWeightPt[iCut] = new TH1F("MC_Meson_WOWeights_Pt", "MC_Meson_WOWeights_Pt; p_{T} (GeV/c)",
1113  (Int_t)((maxPt-minPt)/binWidthPt), minPt, maxPt);
1114  fMCList[iCut]->Add(fHistoMCMesonWOWeightPt[iCut]);
1115 
1116  fHistoMCMesonInAccPt[iCut] = new TH1F("MC_MesonInAcc_Pt", "MC_MesonInAcc_Pt; p_{T} (GeV/c)",
1117  (Int_t)((maxPt-minPt)/binWidthPt), minPt, maxPt);
1118  fHistoMCMesonInAccPt[iCut]->Sumw2();
1119  fMCList[iCut]->Add(fHistoMCMesonInAccPt[iCut]);
1120  fHistoMCMesonWOWeightInAccPt[iCut] = new TH1F("MC_MesonWOWeightInAcc_Pt", "MC_MesonWOWeightInAcc_Pt; p_{T} (GeV/c)",
1121  (Int_t)((maxPt-minPt)/binWidthPt), minPt, maxPt);
1122  fMCList[iCut]->Add(fHistoMCMesonWOWeightInAccPt[iCut]);
1123 
1124  if (fIsMC > 1){
1125  fHistoMCMesonWOWeightPt[iCut]->Sumw2();
1126  fHistoMCMesonWOWeightInAccPt[iCut]->Sumw2();
1127  fHistoMCMesonWOEvtWeightPt[iCut] = new TH1F("MC_Meson_WOEventWeights_Pt", "MC_Meson_WOEventWeights_Pt; p_{T} (GeV/c)",
1128  (Int_t)((maxPt-minPt)/binWidthPt), minPt, maxPt);
1129  fMCList[iCut]->Add(fHistoMCMesonWOEvtWeightPt[iCut]);
1130  fHistoMCMesonWOEvtWeightInAccPt[iCut] = new TH1F("MC_Meson_WOEventWeightsInAcc_Pt", "MC_Meson_WOEventWeightsInAcc_Pt; p_{T} (GeV/c)",
1131  (Int_t)((maxPt-minPt)/binWidthPt), minPt, maxPt);
1132  fMCList[iCut]->Add(fHistoMCMesonWOEvtWeightInAccPt[iCut]);
1133 
1134  if (fDoMesonQA > 0 && fIsMC == 2){
1135  fHistoMCMesonPtJetPt[iCut] = new TH2F("MC_Meson_Pt_JetPt", "MC_Meson_Pt_JetPt; p_{T} (GeV/c); p_{T,jet} (GeV/c)",
1136  nBinsQAPt, arrQAPtBinning, 200, 0, 200);
1137  fHistoMCMesonPtJetPt[iCut]->Sumw2();
1138  fMCList[iCut]->Add(fHistoMCMesonPtJetPt[iCut]);
1139  }
1140  }
1141 
1142  // book additional MC QA histograms
1143  if (fDoMesonQA > 0){
1144  fHistoMCMesonPtY[iCut] = new TH2F("MC_Meson_Pt_Y", "MC_Meson_Pt_Y; p_{T} (GeV/c); Y ",
1145  nBinsQAPt, arrQAPtBinning, 150, -1.5, 1.5);
1146  fHistoMCMesonPtY[iCut]->Sumw2();
1147  fMCList[iCut]->Add(fHistoMCMesonPtY[iCut]);
1148  fHistoMCMesonPtAlpha[iCut] = new TH2F("MC_Meson_Pt_Alpha", "MC_Meson_Pt_Alpha; p_{T} (GeV/c); #alpha",
1149  nBinsQAPt, arrQAPtBinning, 200, -1, 1);
1150  fMCList[iCut]->Add(fHistoMCMesonPtAlpha[iCut]);
1151 
1152  if (fIsMC == 2){
1153  fHistoMCMesonPtAlpha[iCut]->Sumw2();
1154  }
1155  }
1156 
1157  fTrueList[iCut] = new TList();
1158  fTrueList[iCut]->SetName(Form("%s True histograms",fullCutString.Data()));
1159  fTrueList[iCut]->SetOwner(kTRUE);
1160  fCutFolder[iCut]->Add(fTrueList[iCut]);
1161 
1162  if(!fDoLightOutput){
1163  if (fMesonRecoMode < 2){
1164  fHistoTrueConvGammaPt[iCut] = new TH1F("ESD_TrueConvGamma_Pt", "ESD_TrueConvGamma_Pt; p_{T,clus} (GeV/c); counts",
1165  (Int_t)((maxPt-minPt)/binWidthPt), minPt, maxPt);
1166  fTrueList[iCut]->Add(fHistoTrueConvGammaPt[iCut]);
1167  fHistoDoubleCountTrueConvGammaRPt[iCut] = new TH2F("ESD_TrueDoubleCountConvGamma_R_Pt", "ESD_TrueDoubleCountConvGamma_R_Pt; R (cm); p_{T,clus} (GeV/c)",
1168  800, 0, 200, (Int_t)((maxPt-minPt)/binWidthPt), minPt, maxPt);
1169  fTrueList[iCut]->Add(fHistoDoubleCountTrueConvGammaRPt[iCut]);
1170  fHistoMultipleCountTrueConvGamma[iCut] = new TH1F("ESD_TrueMultipleCountConvGamma", "ESD_TrueMultipleCountConvGamma", 10, 1, 11);
1171  fTrueList[iCut]->Add(fHistoMultipleCountTrueConvGamma[iCut]);
1172  fHistoTruePrimaryConvGammaPt[iCut] = new TH1F("ESD_TruePrimaryConvGamma_Pt", "ESD_TruePrimaryConvGamma_Pt;p_{T,clus} (GeV/c); counts", (Int_t)((maxPt-minPt)/binWidthPt), minPt, maxPt);
1173  fTrueList[iCut]->Add(fHistoTruePrimaryConvGammaPt[iCut]);
1174  }
1175  if (fMesonRecoMode > 0){
1176  fHistoTrueClusGammaPt[iCut] = new TH1F("TrueClusGamma_Pt", "ESD_TrueClusGamma_Pt; p_{T,clus} (GeV/c); counts", nBinsClusterPt, arrClusPtBinning);
1177  fTrueList[iCut]->Add(fHistoTrueClusGammaPt[iCut]);
1178  fHistoTruePrimaryClusGammaPt[iCut] = new TH1F("TruePrimaryClusGamma_Pt", "ESD_TruePrimaryClusGamma_Pt; p_{T,clus} (GeV/c); counts", nBinsClusterPt, arrClusPtBinning);
1179  fTrueList[iCut]->Add(fHistoTruePrimaryClusGammaPt[iCut]);
1180  fHistoTrueClusConvGammaPt[iCut] = new TH1F("TrueClusConvGamma_Pt", "TrueClusConvGamma_Pt; p_{T,clus} (GeV/c); counts", nBinsClusterPt, arrClusPtBinning);
1181  fTrueList[iCut]->Add(fHistoTrueClusConvGammaPt[iCut]);
1182  fHistoTruePrimaryClusConvGammaPt[iCut] = new TH1F("TruePrimaryClusConvGamma_Pt", "ESD_TruePrimaryClusConvGamma_Pt; p_{T,clus} (GeV/c); counts", nBinsClusterPt, arrClusPtBinning);
1183  fTrueList[iCut]->Add(fHistoTruePrimaryClusConvGammaPt[iCut]);
1184  fHistoDoubleCountTrueClusterGammaPt[iCut] = new TH2F("TrueDoubleCountClusterGamma_Pt", "TrueDoubleCountClusterGamma_Pt; p_{T,clus} (GeV/c); counts",
1185  nBinsClusterPt, arrClusPtBinning, 2, 0, 2);
1186  fTrueList[iCut]->Add(fHistoDoubleCountTrueClusterGammaPt[iCut]);
1187  fHistoMultipleCountTrueClusterGamma[iCut] = new TH1F("TrueMultipleCountClusterGamma", "TrueMultipleCountClusterGamma", 10, 1, 11);
1188  fTrueList[iCut]->Add(fHistoMultipleCountTrueClusterGamma[iCut]);
1189  fHistoTrueClusConvGammaFullyPt[iCut] = new TH1F("TrueClusConvGammaFullyContained_Pt", "TrueClusConvGammaFullyContained_Pt; p_{T,clus} (GeV/c); counts",
1190  nBinsClusterPt, arrClusPtBinning);
1191  fTrueList[iCut]->Add(fHistoTrueClusConvGammaFullyPt[iCut]);
1192  fHistoTrueNLabelsInClusPt[iCut] = new TH2F("TrueNLabelsInClus_Pt", "TrueNLabelsInClus_Pt; p_{T,clus} (GeV/c); counts",
1193  100, -0.5, 99.5, nBinsClusterPt, arrClusPtBinning);
1194  fTrueList[iCut]->Add(fHistoTrueNLabelsInClusPt[iCut]);
1195  }
1196  }
1197 
1198  if (fIsMC > 1){
1199  if(!fDoLightOutput){
1200  if (fMesonRecoMode < 2){
1201  fHistoTrueConvGammaPt[iCut]->Sumw2();
1202  fHistoDoubleCountTrueConvGammaRPt[iCut]->Sumw2();
1203  fHistoMultipleCountTrueConvGamma[iCut]->Sumw2();
1204  }
1205  if (fMesonRecoMode > 0){
1206  fHistoTruePrimaryConvGammaPt[iCut]->Sumw2();
1207  fHistoTruePrimaryClusGammaPt[iCut]->Sumw2();
1208  fHistoTrueNLabelsInClusPt[iCut]->Sumw2();
1209  fHistoTrueClusConvGammaPt[iCut]->Sumw2();
1210  fHistoTruePrimaryClusConvGammaPt[iCut]->Sumw2();
1211  fHistoTrueClusGammaPt[iCut]->Sumw2();
1212  fHistoDoubleCountTrueClusterGammaPt[iCut]->Sumw2();
1213  fHistoMultipleCountTrueClusterGamma[iCut]->Sumw2();
1214  fHistoTrueClusConvGammaFullyPt[iCut]->Sumw2();
1215  }
1216  }
1217 
1218  }
1219 
1220  fHistoTrueMesonInvMassPt[iCut] = new TH2F("ESD_TrueMeson_InvMass_Pt", "ESD_TrueMeson_InvMass_Pt;M_{inv}(GeV/c^{2});p_{T}(GeV/c)",
1221  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
1222  fHistoTrueMesonInvMassPt[iCut]->Sumw2();
1223  fTrueList[iCut]->Add(fHistoTrueMesonInvMassPt[iCut]);
1224 
1225  fHistoDoubleCountTrueMesonInvMassPt[iCut] = new TH2F("ESD_TrueDoubleCountMeson_InvMass_Pt", "ESD_TrueDoubleCountMeson_InvMass_Pt;M_{inv}(GeV/c^{2});p_{T}(GeV/c)",
1226  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
1227  fTrueList[iCut]->Add(fHistoDoubleCountTrueMesonInvMassPt[iCut]);
1228  fHistoMultipleCountTrueMeson[iCut] = new TH1F("ESD_TrueMultipleCountMeson", "ESD_TrueMultipleCountMeson;# number of multiple counts;#", 10, 1, 11);
1229  fHistoMultipleCountTrueMeson[iCut]->Sumw2();
1230  fTrueList[iCut]->Add(fHistoMultipleCountTrueMeson[iCut]);
1231 
1232  fHistoTruePrimaryMesonInvMassPt[iCut] = new TH2F("ESD_TruePrimaryMeson_InvMass_Pt", "ESD_TruePrimaryMeson_InvMass_Pt;M_{inv}(GeV/c^{2});p_{T}(GeV/c)",
1233  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
1234  fHistoTruePrimaryMesonInvMassPt[iCut]->Sumw2();
1235  fTrueList[iCut]->Add(fHistoTruePrimaryMesonInvMassPt[iCut]);
1236 
1237  fHistoTruePrimaryMesonW0WeightingInvMassPt[iCut] = new TH2F("ESD_TruePrimaryMesonW0Weights_InvMass_Pt",
1238  "ESD_TruePrimaryMesonW0Weights_InvMass_Pt;M_{inv}(GeV/c^{2});p_{T}(GeV/c)",
1239  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
1240  fTrueList[iCut]->Add(fHistoTruePrimaryMesonW0WeightingInvMassPt[iCut]);
1241 
1242  fProfileTruePrimaryMesonWeightsInvMassPt[iCut] = new TProfile2D("ESD_TruePrimaryMesonWeights_InvMass_Pt",
1243  "ESD_TruePrimaryMesonWeights_InvMass_Pt;M_{inv}(GeV/c^{2});p_{T}(GeV/c)",
1244  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
1245  fProfileTruePrimaryMesonWeightsInvMassPt[iCut]->Sumw2();
1246  fTrueList[iCut]->Add(fProfileTruePrimaryMesonWeightsInvMassPt[iCut]);
1247 
1248  if (fMesonRecoMode == 1){
1249  fHistoTrueMesonMatchedInvMassPt[iCut] = new TH2F("ESD_TrueMeson_Matched_InvMass_Pt", "ESD_TrueMeson_Matched_InvMass_Pt;M_{inv}(GeV/c^{2});p_{T}(GeV/c)",
1250  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
1251  fHistoTrueMesonMatchedInvMassPt[iCut]->Sumw2();
1252  fTrueList[iCut]->Add(fHistoTrueMesonMatchedInvMassPt[iCut]);
1253  }
1254 
1255  if (fIsMC > 1){
1256  fHistoTruePrimaryMesonW0WeightingInvMassPt[iCut]->Sumw2();
1257  }
1258 
1259  if (fDoMesonQA > 0){
1260  if (fIsMC < 2){
1261  if (fMesonRecoMode > 0){
1262  fHistoTrueMesonCaloPhotonInvMassPt[iCut] = new TH2F( "ESD_TrueMesonCaloPhoton_InvMass_Pt",
1263  "ESD_TrueMesonCaloPhoton_InvMass_Pt; M_{inv}(GeV/c^{2}) #gamma #gamma;p_{T}(GeV/c) ",
1264  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
1265  fTrueList[iCut]->Add(fHistoTrueMesonCaloPhotonInvMassPt[iCut]);
1266  fHistoTrueMesonCaloConvertedPhotonInvMassPt[iCut] = new TH2F( "ESD_TrueMesonCaloConvertedPhoton_InvMass_Pt",
1267  "ESD_TrueMesonCaloConvertedPhoton_InvMass_Pt;M_{inv}(GeV/c^{2}) #gamma #gamma_{conv};p_{T}(GeV/c)",
1268  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
1269  fTrueList[iCut]->Add(fHistoTrueMesonCaloConvertedPhotonInvMassPt[iCut]);
1270  fHistoTrueMesonCaloElectronInvMassPt[iCut] = new TH2F("ESD_TrueMesonCaloElectron_InvMass_Pt",
1271  "ESD_TrueMesonCaloElectron_InvMass_Pt; M_{inv}(GeV/c^{2}) #gamma e^{#pm}; ; p_{T}(GeV/c)",
1272  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
1273  fTrueList[iCut]->Add(fHistoTrueMesonCaloElectronInvMassPt[iCut]);
1274  fHistoTrueMesonCaloMergedClusterInvMassPt[iCut] = new TH2F("ESD_TrueMesonCaloMergedCluster_InvMass_Pt",
1275  "ESD_TrueMesonCaloMergedCluster_InvMass_Pt; M_{inv}(GeV/c^{2}) #gamma merged cluster; p_{T}(GeV/c)",
1276  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
1277  fTrueList[iCut]->Add(fHistoTrueMesonCaloMergedClusterInvMassPt[iCut]);
1278  fHistoTrueMesonCaloMergedClusterPartConvInvMassPt[iCut] = new TH2F( "ESD_TrueMesonCaloMergedClusterPartConv_InvMass_Pt",
1279  "ESD_TrueMesonCaloMergedClusterPartConv_InvMass_Pt; M_{inv}(GeV/c^{2}) #gamma merged cluster, part conv; p_{T}(GeV/c)",
1280  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
1281  fTrueList[iCut]->Add(fHistoTrueMesonCaloMergedClusterPartConvInvMassPt[iCut]);
1282  }
1283  if (fMesonRecoMode == 1){
1284  fHistoTrueMesonCaloConvertedPhotonMatchedInvMassPt[iCut] = new TH2F("ESD_TrueMesonCaloConvertedPhotonMatched_InvMass_Pt",
1285  "ESD_TrueMesonCaloConvertedPhotonMatched_InvMass_Pt;M_{inv}(GeV/c^{2}) #gamma #gamma_{conv,matched};p_{T}(GeV/c)",
1286  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
1287  fTrueList[iCut]->Add(fHistoTrueMesonCaloConvertedPhotonMatchedInvMassPt[iCut]);
1288  fHistoTrueMotherMesonConvPhotonEtaPhi[iCut] = new TH2F("ESD_TrueMotherMesonConvPhoton_Eta_Phi", "conv photons for true ; #phi_{#gamma_{conv}}(rad);#eta_{#gamma_{conv}}",
1289  600, 0, 2*TMath::Pi(), 200, -1, 1);
1290  fTrueList[iCut]->Add(fHistoTrueMotherMesonConvPhotonEtaPhi[iCut]);
1291  }
1292  if (fMesonRecoMode == 2){
1293  fHistoTrueMesonCaloMixedPhotonConvPhotonInvMassPt[iCut] = new TH2F("ESD_TrueMesonCaloMixedPhotonConvertedPhoton_InvMass_Pt",
1294  "ESD_TrueMesonCaloMixedPhotonConvertedPhoton_InvMass_Pt;M_{inv}(GeV/c^{2}) #gamma #gamma_{conv,matched};p_{T}(GeV/c)",
1295  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
1296  fTrueList[iCut]->Add(fHistoTrueMesonCaloMixedPhotonConvPhotonInvMassPt[iCut]);
1297  }
1298  fHistoTruePrimaryMesonMCPtResolPt[iCut] = new TH2F("ESD_TruePrimaryMeson_MCPt_ResolPt",
1299  "ESD_TruePrimaryMeson_ResolPt_MCPt; p_{T,MC}(GeV/c); (p_{T,rec}-p_{T,MC})/p_{T,MC}()",
1300  500, 0.03, 25, 1000, -1., 1.);
1301  fHistoTruePrimaryMesonMCPtResolPt[iCut]->Sumw2();
1302  SetLogBinningXTH2(fHistoTruePrimaryMesonMCPtResolPt[iCut]);
1303  fTrueList[iCut]->Add(fHistoTruePrimaryMesonMCPtResolPt[iCut]);
1304  }
1305  if(fDoMesonQA>1){
1306  fHistoTrueBckGGInvMassPt[iCut] = new TH2F("ESD_TrueBckGG_InvMass_Pt",
1307  "ESD_TrueBckGG_InvMass_Pt; M_{inv} (GeV/c^{2}) #gamma #gamma no signal; #pair p_{T}(GeV/c)",
1308  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
1309  fTrueList[iCut]->Add(fHistoTrueBckGGInvMassPt[iCut]);
1310  fHistoTrueBckFullMesonContainedInOneClusterInvMassPt[iCut] = new TH2F( "ESD_TrueBckFullMesonContained_InvMass_Pt",
1311  "ESD_TrueBckFullMesonContained_InvMass_Pt; M_{inv} (GeV/c^{2}) #gamma #gamma, calo gamma with full pi0; #pair p_{T}(GeV/c)",
1312  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
1313  fTrueList[iCut]->Add(fHistoTrueBckFullMesonContainedInOneClusterInvMassPt[iCut]);
1314  fHistoTrueBckAsymEClustersInvMassPt[iCut] = new TH2F("ESD_TrueBckAsymEClus_InvMass_Pt",
1315  "ESD_TrueBckAsymEClus_InvMass_Pt; M_{inv} (GeV/c^{2}) #gamma #gamma, calo gamma >70% of pi0 energy; #pair p_{T}(GeV/c)",
1316  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
1317  fTrueList[iCut]->Add(fHistoTrueBckAsymEClustersInvMassPt[iCut]);
1318  fHistoTrueBckContInvMassPt[iCut] = new TH2F("ESD_TrueBckCont_InvMass_Pt",
1319  "ESD_TrueBckCont_InvMass_Pt; M_{inv} (GeV/c^{2}) contamination; #pair p_{T}(GeV/c)",
1320  fMesonInvMassNBins, fMesonInvMassMin, fMesonInvMassMax, nBinsPt, arrPtBinning);
1321  fTrueList[iCut]->Add(fHistoTrueBckContInvMassPt[iCut]);
1322  }
1323  fHistoTrueMesonPtY[iCut] = new TH2F("ESD_TrueMeson_Pt_Y", "ESD_TrueMeson_Pt_Y; p_{T}(GeV/c); Y",
1324  nBinsQAPt, arrQAPtBinning, 150, -1.5, 1.5);
1325  fTrueList[iCut]->Add(fHistoTrueMesonPtY[iCut]);
1326  fHistoTrueMesonPtAlpha[iCut] = new TH2F("ESD_TrueMeson_Pt_Alpha", "ESD_TrueMeson_Pt_Alpha; p_{T}(GeV/c); #alpha",
1327  nBinsQAPt, arrQAPtBinning, 200, -1, 1);
1328  fTrueList[iCut]->Add(fHistoTrueMesonPtAlpha[iCut]);
1329  fHistoTrueMesonPtOpenAngle[iCut] = new TH2F("ESD_TrueMeson_Pt_OpenAngle", "ESD_TrueMeson_Pt_OpenAngle; p_{T}(GeV/c); #theta",
1330  nBinsQAPt, arrQAPtBinning, 100, 0, 1);
1331  fTrueList[iCut]->Add(fHistoTrueMesonPtOpenAngle[iCut]);
1332 
1333  if (fIsMC > 1){
1334  if(fDoMesonQA>1){
1335  fHistoTrueBckGGInvMassPt[iCut]->Sumw2();
1336  fHistoTrueBckFullMesonContainedInOneClusterInvMassPt[iCut]->Sumw2();
1337  fHistoTrueBckAsymEClustersInvMassPt[iCut]->Sumw2();
1338  fHistoTrueBckContInvMassPt[iCut]->Sumw2();
1339  }
1340  fHistoTrueMesonPtY[iCut]->Sumw2();
1341  fHistoTrueMesonPtAlpha[iCut]->Sumw2();
1342  fHistoTrueMesonPtOpenAngle[iCut]->Sumw2();
1343  if (fHistoTrueMotherMesonConvPhotonEtaPhi[iCut]) fHistoTrueMotherMesonConvPhotonEtaPhi[iCut]->Sumw2();
1344  }
1345  }
1346  }
1347  }
1348 
1349  fVectorDoubleCountTrueMesons.clear();
1350  fVectorDoubleCountTrueConvGammas.clear();
1351  fVectorDoubleCountTrueClusterGammas.clear();
1352 
1353  fMapMultipleCountTrueMesons.clear();
1354  fMapMultipleCountTrueConvGammas.clear();
1355  fMapMultipleCountTrueClusterGammas.clear();
1356 
1357  fVectorRecTrueMesons.clear();
1358 
1359  if(fV0Reader)
1360  if((AliConversionPhotonCuts*)fV0Reader->GetConversionCuts())
1361  if(((AliConversionPhotonCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms())
1362  fOutputContainer->Add(((AliConversionPhotonCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms());
1363  if(fV0Reader)
1364  if((AliConvEventCuts*)fV0Reader->GetEventCuts())
1365  if(((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetCutHistograms())
1366  fOutputContainer->Add(((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetCutHistograms());
1367 
1368  if( fV0Reader && fV0Reader->GetProduceV0FindingEfficiency())
1369  if (fV0Reader->GetV0FindingEfficiencyHistograms())
1370  fOutputContainer->Add(fV0Reader->GetV0FindingEfficiencyHistograms());
1371 
1372  for(Int_t iMatcherTask = 0; iMatcherTask < 3; iMatcherTask++){
1373  AliCaloTrackMatcher* temp = (AliCaloTrackMatcher*) (AliAnalysisManager::GetAnalysisManager()->GetTask(Form("CaloTrackMatcher_%i",iMatcherTask)));
1374  if(temp) fOutputContainer->Add(temp->GetCaloTrackMatcherHistograms());
1375  }
1376 
1377  for(Int_t iCut = 0; iCut<fnCuts;iCut++){
1378  if(!((AliConvEventCuts*)fEventCutArray->At(iCut))) continue;
1379  if(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutHistograms()){
1380  fCutFolder[iCut]->Add(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutHistograms());
1381  }
1382  if (fMesonRecoMode < 2){
1383  if(!((AliConversionPhotonCuts*)fCutArray->At(iCut))) continue;
1384  if(((AliConversionPhotonCuts*)fCutArray->At(iCut))->GetCutHistograms()){
1385  fCutFolder[iCut]->Add(((AliConversionPhotonCuts*)fCutArray->At(iCut))->GetCutHistograms());
1386  }
1387  }
1388  if (fMesonRecoMode > 0){
1389  if(!((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))) continue;
1390  if(((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutHistograms()){
1391  fCutFolder[iCut]->Add(((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutHistograms());
1392  }
1393  if(fSetPlotHistsExtQA){
1394  if(!((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))) continue;
1395  if(((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetExtQAHistograms()){
1396  fCutFolder[iCut]->Add(((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetExtQAHistograms());
1397  }
1398  }
1399  }
1400  if(!((AliConversionMesonCuts*)fMesonCutArray->At(iCut))) continue;
1401  if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutHistograms()){
1402  fCutFolder[iCut]->Add(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutHistograms());
1403  }
1404  }
1405 
1406  if (fIsMC > 0){
1407  tBrokenFiles = new TTree("BrokenFiles", "BrokenFiles");
1408  tBrokenFiles->Branch("fileName",&fFileNameBroken);
1409  fOutputContainer->Add(tBrokenFiles);
1410  }
1411 
1412  PostData(1, fOutputContainer);
1413 }
1414 //_____________________________________________________________________________
1415 Bool_t AliAnalysisTaskHeavyNeutralMesonToGG::Notify(){
1416 
1417  for(Int_t iCut = 0; iCut<fnCuts;iCut++){
1418  if (((AliConvEventCuts*)fEventCutArray->At(iCut))->GetPeriodEnum() == AliConvEventCuts::kNoPeriod && ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetPeriodEnum() != AliConvEventCuts::kNoPeriod){
1419  ((AliConvEventCuts*)fEventCutArray->At(iCut))->SetPeriodEnumExplicit(((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetPeriodEnum());
1420  } else if (((AliConvEventCuts*)fEventCutArray->At(iCut))->GetPeriodEnum() == AliConvEventCuts::kNoPeriod ){
1421  ((AliConvEventCuts*)fEventCutArray->At(iCut))->SetPeriodEnum(fV0Reader->GetPeriodName());
1422  }
1423 
1424  if(fIsHeavyIon == 2){
1425  if(!((AliConvEventCuts*)fEventCutArray->At(iCut))->GetDoEtaShift()){
1426  fProfileEtaShift[iCut]->Fill(0.,(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift()));
1427  continue; // No Eta Shift requested, continue
1428  }
1429  if(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift() == 0.0){ // Eta Shift requested but not set, get shift automatically
1430  ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCorrectEtaShiftFromPeriod();
1431  fProfileEtaShift[iCut]->Fill(0.,(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift()));
1432  ((AliConvEventCuts*)fEventCutArray->At(iCut))->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
1433  continue;
1434  } else{
1435  printf(" Gamma Conversion Task %s :: Eta Shift Manually Set to %f \n\n",
1436  (((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber()).Data(),((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift());
1437  fProfileEtaShift[iCut]->Fill(0.,(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift()));
1438  ((AliConvEventCuts*)fEventCutArray->At(iCut))->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
1439  }
1440  }
1441  }
1442 
1443  return kTRUE;
1444 }
1445 //_____________________________________________________________________________
1446 void AliAnalysisTaskHeavyNeutralMesonToGG::UserExec(Option_t *){
1447  //
1448  // Called for each event
1449  //
1450  fInputEvent = InputEvent();
1451 
1452  if(fIsMC > 0) fMCEvent = MCEvent();
1453 
1454  Int_t eventQuality = ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEventQuality();
1455  if(fInputEvent->IsIncompleteDAQ()==kTRUE) eventQuality = 2;// incomplete event
1456  // Event Not Accepted due to MC event missing or wrong trigger for V0ReaderV1 or because it is incomplete => abort processing of this event/file
1457  if(eventQuality == 2 || eventQuality == 3){
1458  // write out name of broken file for first event
1459  if (fIsMC > 0){
1460  if (fInputEvent->IsA()==AliESDEvent::Class()){
1461  if (((AliESDEvent*)fInputEvent)->GetEventNumberInFile() == 0){
1462  fFileNameBroken = new TObjString(Form("%s",((TString)fV0Reader->GetCurrentFileName()).Data()));
1463  if (tBrokenFiles) tBrokenFiles->Fill();
1464  delete fFileNameBroken;
1465  }
1466  }
1467  }
1468 
1469  for(Int_t iCut = 0; iCut<fnCuts; iCut++){
1470  fHistoNEvents[iCut]->Fill(eventQuality);
1471  if (fIsMC > 1) fHistoNEventsWOWeight[iCut]->Fill(eventQuality);
1472  }
1473  return;
1474  }
1475 
1476  if(fInputEvent->IsA()==AliAODEvent::Class()){
1477  fInputEvent->InitMagneticField();
1478  }
1479 
1480  fReaderGammas = fV0Reader->GetReconstructedGammas(); // Gammas from default Cut
1481 
1482  // ------------------- BeginEvent ----------------------------
1483  AliEventplane *EventPlane = fInputEvent->GetEventplane();
1484  if(fIsHeavyIon ==1)fEventPlaneAngle = EventPlane->GetEventplane("V0",fInputEvent,2);
1485  else fEventPlaneAngle=0.0;
1486 
1487  if(fIsMC>0 && fInputEvent->IsA()==AliAODEvent::Class() && !(fV0Reader->AreAODsRelabeled()) && fMesonRecoMode < 2){
1488  RelabelAODPhotonCandidates(kTRUE);// In case of AODMC relabeling MC
1489  fV0Reader->RelabelAODs(kTRUE);
1490  }
1491 
1492  for(Int_t iCut = 0; iCut<fnCuts; iCut++){
1493 
1494  fiCut = iCut;
1495  Bool_t isRunningEMCALrelAna = kFALSE;
1496  if (fMesonRecoMode > 0){
1497  if (((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->GetClusterType() == 1) isRunningEMCALrelAna = kTRUE;
1498  }
1499 
1500  Int_t eventNotAccepted = ((AliConvEventCuts*)fEventCutArray->At(iCut))->IsEventAcceptedByCut(fV0Reader->GetEventCuts(),fInputEvent,fMCEvent,fIsHeavyIon,isRunningEMCALrelAna);
1501 
1502  if(fIsMC==2){
1503  Float_t xsection = -1.;
1504  Float_t ntrials = -1.;
1505  ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetXSectionAndNTrials(fMCEvent,xsection,ntrials,fInputEvent);
1506  if((xsection==-1.) || (ntrials==-1.)) AliFatal(Form("ERROR: GetXSectionAndNTrials returned invalid xsection/ntrials, periodName from V0Reader: '%s'",fV0Reader->GetPeriodName().Data()));
1507  fProfileJetJetXSection[iCut]->Fill(0.,xsection);
1508  fHistoJetJetNTrials[iCut]->Fill("#sum{NTrials}",ntrials);
1509  }
1510 
1511  if(fIsMC>0){
1512  fWeightJetJetMC = 1;
1513  // cout << fMCEvent << endl;
1514  Bool_t isMCJet = ((AliConvEventCuts*)fEventCutArray->At(iCut))->IsJetJetMCEventAccepted( fMCEvent, fWeightJetJetMC, fInputEvent );
1515  if (fIsMC == 3){
1516  Double_t weightMult = ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetWeightForMultiplicity(fV0Reader->GetNumberOfPrimaryTracks());
1517  fWeightJetJetMC = fWeightJetJetMC*weightMult;
1518  }
1519  if(fIsMC==1) fWeightJetJetMC = 1;
1520  if (!isMCJet){
1521  fHistoNEvents[iCut]->Fill(10,fWeightJetJetMC);
1522  if (fIsMC>1) fHistoNEventsWOWeight[iCut]->Fill(10);
1523  continue;
1524  }
1525  }
1526 
1527  Bool_t triggered = kTRUE;
1528 
1529  if(eventNotAccepted!= 0){
1530  // cout << "event rejected due to wrong trigger: " <<eventNotAccepted << endl;
1531  fHistoNEvents[iCut]->Fill(eventNotAccepted, fWeightJetJetMC); // Check Centrality, PileUp, SDD and V0AND --> Not Accepted => eventQuality = 1
1532  if (fIsMC>1) fHistoNEventsWOWeight[iCut]->Fill(eventNotAccepted);
1533  if (eventNotAccepted==3 && fIsMC > 0){
1534  triggered = kFALSE;
1535  }else {
1536  continue;
1537  }
1538  }
1539 
1540  if(eventQuality != 0){// Event Not Accepted
1541  //cout << "event rejected due to: " <<eventQuality << endl;
1542  fHistoNEvents[iCut]->Fill(eventQuality, fWeightJetJetMC);
1543  if (fIsMC>1) fHistoNEventsWOWeight[iCut]->Fill(eventQuality);
1544  continue;
1545  }
1546 
1547  if (triggered==kTRUE){
1548  fHistoNEvents[iCut]->Fill(eventQuality, fWeightJetJetMC); // Should be 0 here
1549  if (fIsMC>1) fHistoNEventsWOWeight[iCut]->Fill(eventQuality); // Should be 0 here
1550 
1551  fHistoNGoodESDTracks[iCut]->Fill(fV0Reader->GetNumberOfPrimaryTracks(), fWeightJetJetMC);
1552  fHistoVertexZ[iCut]->Fill(fInputEvent->GetPrimaryVertex()->GetZ(), fWeightJetJetMC);
1553  if(!fDoLightOutput){
1554  fHistoVertexX[iCut]->Fill(fInputEvent->GetPrimaryVertex()->GetX(), fWeightJetJetMC);
1555  fHistoVertexY[iCut]->Fill(fInputEvent->GetPrimaryVertex()->GetY(), fWeightJetJetMC);
1556  fHistoSPDClusterTrackletBackground[iCut]->Fill(fInputEvent->GetMultiplicity()->GetNumberOfTracklets(),(fInputEvent->GetNumberOfITSClusters(0)+fInputEvent->GetNumberOfITSClusters(1)),fWeightJetJetMC);
1557  if(((AliConvEventCuts*)fEventCutArray->At(iCut))->IsHeavyIon() == 2) fHistoNV0Tracks[iCut]->Fill(fInputEvent->GetVZEROData()->GetMTotV0A(), fWeightJetJetMC);
1558  else fHistoNV0Tracks[iCut]->Fill(fInputEvent->GetVZEROData()->GetMTotV0A()+fInputEvent->GetVZEROData()->GetMTotV0C(), fWeightJetJetMC);
1559  }
1560  }
1561 
1562  if(fIsMC>0){
1563  // Process MC Particle
1564  if(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetSignalRejection() != 0){
1565  if(fInputEvent->IsA()==AliESDEvent::Class()){
1566  ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetNotRejectedParticles(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetSignalRejection(),
1567  ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetAcceptedHeader(),
1568  fMCEvent);
1569  }
1570  else if(fInputEvent->IsA()==AliAODEvent::Class()){
1571  ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetNotRejectedParticles(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetSignalRejection(),
1572  ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetAcceptedHeader(),
1573  fInputEvent);
1574  }
1575 
1576  if(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetAcceptedHeader()){
1577  for(Int_t i = 0;i<(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetAcceptedHeader())->GetEntries();i++){
1578  TString nameBin= fHistoMCHeaders[iCut]->GetXaxis()->GetBinLabel(i+1);
1579  if (nameBin.CompareTo("")== 0){
1580  TString nameHeader = ((TObjString*)((TList*)((AliConvEventCuts*)fEventCutArray->At(iCut))
1581  ->GetAcceptedHeader())->At(i))->GetString();
1582  fHistoMCHeaders[iCut]->GetXaxis()->SetBinLabel(i+1,nameHeader.Data());
1583  }
1584  }
1585  }
1586  }
1587  }
1588  if(fIsMC>0){
1589  if(fInputEvent->IsA()==AliESDEvent::Class())
1590  ProcessMCParticles();
1591  if(fInputEvent->IsA()==AliAODEvent::Class())
1592  ProcessAODMCParticles();
1593  }
1594 
1595  if (triggered==kFALSE) continue;
1596 
1597  // it is in the loop to have the same conversion cut string (used also for MC stuff that should be same for V0 and Cluster)
1598  if (fMesonRecoMode > 0) // process calo clusters
1599  ProcessClusters();
1600  if (fMesonRecoMode < 2) // Process this cuts gammas
1601  ProcessPhotonCandidates();
1602 
1603  if (fMesonRecoMode < 2) fHistoNGammaConvCandidates[iCut]->Fill(fGammaCandidates->GetEntries(),fWeightJetJetMC);
1604  if (fMesonRecoMode > 0) fHistoNGammaCaloCandidates[iCut]->Fill(fClusterCandidates->GetEntries(),fWeightJetJetMC);
1605 
1606  if (fMesonRecoMode < 2){
1607  if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseMCPSmearing() && fIsMC>0){
1608  fUnsmearedPx = new Double_t[fGammaCandidates->GetEntries()]; // Store unsmeared Momenta
1609  fUnsmearedPy = new Double_t[fGammaCandidates->GetEntries()];
1610  fUnsmearedPz = new Double_t[fGammaCandidates->GetEntries()];
1611  fUnsmearedE = new Double_t[fGammaCandidates->GetEntries()];
1612 
1613  for(Int_t gamma=0;gamma<fGammaCandidates->GetEntries();gamma++){ // Smear the AODPhotons in MC
1614  fUnsmearedPx[gamma] = ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->Px();
1615  fUnsmearedPy[gamma] = ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->Py();
1616  fUnsmearedPz[gamma] = ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->Pz();
1617  fUnsmearedE[gamma] = ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->E();
1618  ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->SmearParticle(dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(gamma)));
1619  }
1620  }
1621  }
1622 
1623  CalculateMesonCandidates(); // Combine Gammas from conversion and from calo
1624 
1625  if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->DoBGCalculation()){
1626  if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->BackgroundHandlerType() == 0){
1627  CalculateBackground(); // Combinatorial Background
1628  UpdateEventByEventData(); // Store Event for mixed Events
1629  } else{
1630  CalculateBackgroundRP(); // Combinatorial Background
1631  fBGHandlerRP[iCut]->AddEvent(fGammaCandidates,fInputEvent); // Store Event for mixed Events
1632  fBGClusHandlerRP[iCut]->AddEvent(fClusterCandidates,fInputEvent); // Store Event for mixed Events
1633  }
1634  }
1635 
1636  if (fMesonRecoMode < 2){
1637  if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseMCPSmearing() && fIsMC>0){
1638  for(Int_t gamma=0;gamma<fGammaCandidates->GetEntries();gamma++){ // Smear the AODPhotons in MC
1639  ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetPx(fUnsmearedPx[gamma]); // Reset Unsmeared Momenta
1640  ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetPy(fUnsmearedPy[gamma]);
1641  ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetPz(fUnsmearedPz[gamma]);
1642  ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetE(fUnsmearedE[gamma]);
1643  }
1644  delete[] fUnsmearedPx; fUnsmearedPx = 0x0;
1645  delete[] fUnsmearedPy; fUnsmearedPy = 0x0;
1646  delete[] fUnsmearedPz; fUnsmearedPz = 0x0;
1647  delete[] fUnsmearedE;fUnsmearedE = 0x0;
1648  }
1649  }
1650 
1651  if(fIsMC>0){
1652  fVectorRecTrueMesons.clear();
1653  fVectorDoubleCountTrueMesons.clear();
1654  FillMultipleCountHistoAndClear(fMapMultipleCountTrueMesons,fHistoMultipleCountTrueMeson[iCut]);
1655  fVectorDoubleCountTrueConvGammas.clear();
1656  if(!fDoLightOutput && fMesonRecoMode < 2) FillMultipleCountHistoAndClear(fMapMultipleCountTrueConvGammas,fHistoMultipleCountTrueConvGamma[iCut]);
1657  fVectorDoubleCountTrueClusterGammas.clear();
1658  if(!fDoLightOutput && fMesonRecoMode > 0) FillMultipleCountHistoAndClear(fMapMultipleCountTrueClusterGammas,fHistoMultipleCountTrueClusterGamma[iCut]);
1659  }
1660 
1661  fGammaCandidates->Clear(); // delete this cuts good gammas
1662  fClusterCandidates->Clear(); // delete cluster candidates
1663  }
1664 
1665  if(fIsMC>0 && fInputEvent->IsA()==AliAODEvent::Class() && !(fV0Reader->AreAODsRelabeled()) && fMesonRecoMode < 2){
1666  RelabelAODPhotonCandidates(kFALSE); // Back to ESDMC Label
1667  fV0Reader->RelabelAODs(kFALSE);
1668  }
1669 
1670  PostData(1, fOutputContainer);
1671 }
1672 
1673 //________________________________________________________________________
1674 void AliAnalysisTaskHeavyNeutralMesonToGG::ProcessClusters(){
1675  Int_t nclus = 0;
1676  TClonesArray * arrClustersProcess = NULL;
1677  if(!fCorrTaskSetting.CompareTo("")){
1678  nclus = fInputEvent->GetNumberOfCaloClusters();
1679  } else {
1680  arrClustersProcess = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(Form("%sClustersBranch",fCorrTaskSetting.Data())));
1681  if(!arrClustersProcess)
1682  AliFatal(Form("%sClustersBranch was not found in AliAnalysisTaskHeavyNeutralMesonToGG! Check the correction framework settings!",fCorrTaskSetting.Data()));
1683  nclus = arrClustersProcess->GetEntries();
1684  }
1685 
1686  // cout << nclus << endl;
1687 
1688  if(nclus == 0) return;
1689 
1690  // plotting histograms on cell/tower level, only if extendedMatchAndQA > 1
1691  ((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->FillHistogramsExtendedQA(fInputEvent,fIsMC);
1692 
1693  // match tracks to clusters
1694  if(fDoPrimaryTrackMatching) ((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->MatchTracksToClusters(fInputEvent,fWeightJetJetMC,kFALSE);
1695 
1696  // vertex
1697  Double_t vertex[3] = {0,0,0};
1698  InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1699 
1700  // Loop over EMCal clusters
1701  for(Int_t i = 0; i < nclus; i++){
1702  AliVCluster* clus = NULL;
1703  if(fInputEvent->IsA()==AliESDEvent::Class()){
1704  if(arrClustersProcess)
1705  clus = new AliESDCaloCluster(*(AliESDCaloCluster*)arrClustersProcess->At(i));
1706  else
1707  clus = new AliESDCaloCluster(*(AliESDCaloCluster*)fInputEvent->GetCaloCluster(i));
1708  } else if(fInputEvent->IsA()==AliAODEvent::Class()){
1709  if(arrClustersProcess)
1710  clus = new AliAODCaloCluster(*(AliAODCaloCluster*)arrClustersProcess->At(i));
1711  else
1712  clus = new AliAODCaloCluster(*(AliAODCaloCluster*)fInputEvent->GetCaloCluster(i));
1713  }
1714 
1715  if (!clus) continue;
1716  if(!((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->ClusterIsSelected(clus,fInputEvent,fMCEvent,fIsMC,fWeightJetJetMC,i)){
1717  delete clus;
1718  continue;
1719  }
1720 
1721  // TLorentzvector with cluster
1722  TLorentzVector clusterVector;
1723  clus->GetMomentum(clusterVector,vertex);
1724 
1725  TLorentzVector* tmpvec = new TLorentzVector();
1726  tmpvec->SetPxPyPzE(clusterVector.Px(),clusterVector.Py(),clusterVector.Pz(),clusterVector.E());
1727 
1728  // convert to AODConversionPhoton
1729  AliAODConversionPhoton *PhotonCandidate = new AliAODConversionPhoton(tmpvec);
1730  if(!PhotonCandidate){ delete clus; delete tmpvec; continue;}
1731 
1732  // Flag Photon as CaloPhoton
1733  PhotonCandidate->SetIsCaloPhoton();
1734  PhotonCandidate->SetCaloClusterRef((Long_t)i);
1735  // get MC label
1736  if(fIsMC>0){
1737  PhotonCandidate->SetNCaloPhotonMCLabels(clus->GetNLabels());
1738  Int_t* mclabelsCluster = clus->GetLabels();
1739  PhotonCandidate->SetNCaloPhotonMCLabels(clus->GetNLabels());
1740  // cout << clus->GetNLabels() << endl;
1741  if (clus->GetNLabels()>0){
1742  for (Int_t k =0; k<(Int_t)clus->GetNLabels(); k++){
1743  if (k<50)PhotonCandidate->SetCaloPhotonMCLabel(k,mclabelsCluster[k]);
1744  // Int_t pdgCode = fMCEvent->Particle(mclabelsCluster[k])->GetPdgCode();
1745  // cout << "label " << k << "\t" << mclabelsCluster[k] << " pdg code: " << pdgCode << endl;
1746  }
1747  }
1748  }
1749 
1750  fIsFromDesiredHeader = kTRUE;
1751  fIsOverlappingWithOtherHeader = kFALSE;
1752  //TString periodName = fV0Reader->GetPeriodName();
1753  // test whether largest contribution to cluster orginates in added signals
1754  if (fIsMC>0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() > 0){
1755  if (((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetCaloPhotonMCLabel(0), fMCEvent, fInputEvent) == 0){
1756  fIsFromDesiredHeader = kFALSE;
1757  }
1758  if (clus->GetNLabels()>1){
1759  Int_t* mclabelsCluster = clus->GetLabels();
1760  for (Int_t l = 1; l < (Int_t)clus->GetNLabels(); l++ ){
1761  if (((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(mclabelsCluster[l], fMCEvent, fInputEvent) == 0){
1762  fIsOverlappingWithOtherHeader = kTRUE;
1763  }
1764  }
1765  }
1766  }
1767  if (fIsMC > 0 && fHistoClusAllHeadersGammaPt[fiCut])
1768  fHistoClusAllHeadersGammaPt[fiCut]->Fill(PhotonCandidate->Pt(),fWeightJetJetMC);
1769  if (!fIsFromDesiredHeader && fIsMC > 0 && fHistoClusRejectedHeadersGammaPt[fiCut])
1770  fHistoClusRejectedHeadersGammaPt[fiCut]->Fill(PhotonCandidate->Pt(),fWeightJetJetMC);
1771  if (fIsFromDesiredHeader && fIsOverlappingWithOtherHeader && fIsMC > 0 && fHistoClusOverlapHeadersGammaPt[fiCut])
1772  fHistoClusOverlapHeadersGammaPt[fiCut]->Fill(PhotonCandidate->Pt(),fWeightJetJetMC);
1773 
1774  if ( (fIsFromDesiredHeader && !fIsOverlappingWithOtherHeader && !fAllowOverlapHeaders) || (fIsFromDesiredHeader && fAllowOverlapHeaders)){
1775  fHistoClusGammaPt[fiCut]->Fill(PhotonCandidate->Pt(),fWeightJetJetMC);
1776  fHistoClusGammaE[fiCut]->Fill(PhotonCandidate->E(),fWeightJetJetMC);
1777  if(fIsMC>0){
1778  if(fInputEvent->IsA()==AliESDEvent::Class()){
1779  ProcessTrueClusterCandidates(PhotonCandidate,clus->GetM02());
1780  }else {
1781  ProcessTrueClusterCandidatesAOD(PhotonCandidate,clus->GetM02());
1782  }
1783  }
1784  fClusterCandidates->Add(PhotonCandidate); // if no second loop is required add to events good gammas
1785  }else{
1786  delete PhotonCandidate;
1787  }
1788 
1789  delete clus;
1790  delete tmpvec;
1791  }
1792 }
1793 
1794 //________________________________________________________________________
1795 void AliAnalysisTaskHeavyNeutralMesonToGG::ProcessTrueClusterCandidates(AliAODConversionPhoton *TruePhotonCandidate, Float_t clusM02){
1796 
1797  TParticle *Photon = NULL;
1798  if (!TruePhotonCandidate->GetIsCaloPhoton()) AliFatal("CaloPhotonFlag has not been set task will abort");
1799  if (TruePhotonCandidate->GetCaloPhotonMCLabel(0) < 0) return;
1800  if(!fDoLightOutput) fHistoTrueNLabelsInClusPt[fiCut]->Fill(TruePhotonCandidate->GetNCaloPhotonMCLabels(),TruePhotonCandidate->Pt(),fWeightJetJetMC);
1801 
1802  if (TruePhotonCandidate->GetNCaloPhotonMCLabels()>0)Photon = fMCEvent->Particle(TruePhotonCandidate->GetCaloPhotonMCLabel(0));
1803  else return;
1804 
1805  if(Photon == NULL){
1806  // cout << "no photon" << endl;
1807  return;
1808  }
1809 
1810  TruePhotonCandidate->SetCaloPhotonMCFlags(fMCEvent, fEnableSortForClusMC);
1811 
1812  const AliVVertex* primVtxMC = fMCEvent->GetPrimaryVertex();
1813  Double_t mcProdVtxX = primVtxMC->GetX();
1814  Double_t mcProdVtxY = primVtxMC->GetY();
1815  Double_t mcProdVtxZ = primVtxMC->GetZ();
1816  Bool_t isPrimary = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsConversionPrimaryESD( fMCEvent, TruePhotonCandidate->GetCaloPhotonMCLabel(0), mcProdVtxX, mcProdVtxY, mcProdVtxZ);
1817 
1818  // to get primary distrinction right put mother of conversion electron as particle to check
1819  if (TruePhotonCandidate->IsLargestComponentElectron() && TruePhotonCandidate->IsConversion())
1820  isPrimary = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsConversionPrimaryESD( fMCEvent, Photon->GetMother(0), mcProdVtxX, mcProdVtxY, mcProdVtxZ);
1821 
1822 
1823  // Fill histograms for inclusive gamma corrections
1824  // --> a) all clusters with leading real or converted photons
1825  if (TruePhotonCandidate->IsLargestComponentPhoton() || (TruePhotonCandidate->IsLargestComponentElectron() && TruePhotonCandidate->IsConversion()) ){
1826  fHistoTrueClusGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt(),fWeightJetJetMC);
1827  if(!fDoLightOutput){
1828  // how many of those clusters are from converted photons
1829  if (TruePhotonCandidate->IsLargestComponentElectron() && TruePhotonCandidate->IsConversion()){
1830  fHistoTrueClusConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt(), fWeightJetJetMC);
1831  }
1832  // --> b) which of these are primary
1833  if(isPrimary){
1834  fHistoTruePrimaryClusGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt(),fWeightJetJetMC);
1835  // --> c) which are from conversions? Just additonal information
1836  if (TruePhotonCandidate->IsLargestComponentElectron() && TruePhotonCandidate->IsConversion() && (Photon->GetMother(0)>-1)){
1837  fHistoTruePrimaryClusConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt(), fWeightJetJetMC);
1838  }
1839  // --> d) how do the secondaries look like
1840  }
1841  }
1842  }
1843 
1844  // Some additional QA
1845  if (fDoClusterQA > 0){
1846  // how many of the converted photons are fully contained in the cluster
1847  if (TruePhotonCandidate->IsLargestComponentElectron() && TruePhotonCandidate->IsConversion() && TruePhotonCandidate->IsConversionFullyContained())
1848  fHistoTrueClusConvGammaFullyPt[fiCut]->Fill(TruePhotonCandidate->Pt(), fWeightJetJetMC);
1849  }
1850 
1851  // Check if we are double counting photons
1852  Int_t motherLab = Photon->GetMother(0);
1853  if (motherLab > -1){
1854  if (TruePhotonCandidate->IsLargestComponentPhoton()){
1855  if (CheckVectorForDoubleCount(fVectorDoubleCountTrueClusterGammas,motherLab)){
1856  fHistoDoubleCountTrueClusterGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt(),(Double_t)0,fWeightJetJetMC);
1857  FillMultipleCountMap(fMapMultipleCountTrueClusterGammas,motherLab);
1858  }
1859  }
1860  Int_t grandMotherLab = fMCEvent->Particle(motherLab)->GetMother(0);
1861  if (grandMotherLab > -1){
1862  if (TruePhotonCandidate->IsLargestComponentElectron() && TruePhotonCandidate->IsConversion()){
1863  if (CheckVectorForDoubleCount(fVectorDoubleCountTrueClusterGammas,grandMotherLab)){
1864  fHistoDoubleCountTrueClusterGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt(),(Double_t)1,fWeightJetJetMC);
1865  FillMultipleCountMap(fMapMultipleCountTrueClusterGammas,grandMotherLab);
1866  }
1867  }
1868  }
1869  }
1870  return;
1871 }
1872 
1873 
1874 //________________________________________________________________________
1875 void AliAnalysisTaskHeavyNeutralMesonToGG::ProcessTrueClusterCandidatesAOD(AliAODConversionPhoton *TruePhotonCandidate, Float_t clusM02){
1876 
1877  AliAODMCParticle *Photon = NULL;
1878  TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1879  if (AODMCTrackArray){
1880  if (!TruePhotonCandidate->GetIsCaloPhoton()) AliFatal("CaloPhotonFlag has not been set task will abort");
1881  if (TruePhotonCandidate->GetNCaloPhotonMCLabels()>0) Photon = (AliAODMCParticle*) AODMCTrackArray->At(TruePhotonCandidate->GetCaloPhotonMCLabel(0));
1882  else return;
1883  }else {
1884  AliInfo("AODMCTrackArray could not be loaded");
1885  return;
1886  }
1887 
1888  if(Photon == NULL){
1889  // cout << "no photon" << endl;
1890  return;
1891  }
1892 
1893  TruePhotonCandidate->SetCaloPhotonMCFlagsAOD(fInputEvent, fEnableSortForClusMC);
1894  if(!fDoLightOutput) fHistoTrueNLabelsInClusPt[fiCut]->Fill(TruePhotonCandidate->GetNCaloPhotonMCLabels(),TruePhotonCandidate->Pt(),fWeightJetJetMC);
1895 
1896  const AliVVertex* primVtxMC = fMCEvent->GetPrimaryVertex();
1897  Double_t mcProdVtxX = primVtxMC->GetX();
1898  Double_t mcProdVtxY = primVtxMC->GetY();
1899  Double_t mcProdVtxZ = primVtxMC->GetZ();
1900  Bool_t isPrimary = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsConversionPrimaryAOD(fInputEvent, Photon, mcProdVtxX, mcProdVtxY, mcProdVtxZ);
1901 
1902  // to get primary distrinction right put mother of conversion electron as particle to check
1903  if (TruePhotonCandidate->IsLargestComponentElectron() && TruePhotonCandidate->IsConversion()){
1904  if (Photon->GetMother()> -1){
1905  AliAODMCParticle *Mother = (AliAODMCParticle*) AODMCTrackArray->At(Photon->GetMother());
1906  isPrimary = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsConversionPrimaryAOD( fInputEvent, Mother, mcProdVtxX, mcProdVtxY, mcProdVtxZ);
1907  }
1908  }
1909 
1910  // True Photon
1911  if (TruePhotonCandidate->IsLargestComponentPhoton() || (TruePhotonCandidate->IsLargestComponentElectron() && TruePhotonCandidate->IsConversion()) ){
1912  fHistoTrueClusGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt(),fWeightJetJetMC);
1913  if(!fDoLightOutput){
1914  if (TruePhotonCandidate->IsLargestComponentElectron() && TruePhotonCandidate->IsConversion()){
1915  fHistoTrueClusConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt(), fWeightJetJetMC);
1916  }
1917  if(isPrimary){
1918  fHistoTruePrimaryClusGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt(),fWeightJetJetMC);
1919  if (TruePhotonCandidate->IsLargestComponentElectron() && TruePhotonCandidate->IsConversion()){
1920  fHistoTruePrimaryClusConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt(), fWeightJetJetMC);
1921  }
1922  }
1923  }
1924  }
1925  if (fDoClusterQA > 0){
1926  if (TruePhotonCandidate->IsLargestComponentElectron() && TruePhotonCandidate->IsConversion() && TruePhotonCandidate->IsConversionFullyContained())
1927  fHistoTrueClusConvGammaFullyPt[fiCut]->Fill(TruePhotonCandidate->Pt(), fWeightJetJetMC);
1928  }
1929  Int_t motherLab = Photon->GetMother();
1930  if (motherLab > -1){
1931  if (TruePhotonCandidate->IsLargestComponentPhoton()){
1932  if (CheckVectorForDoubleCount(fVectorDoubleCountTrueClusterGammas,motherLab)){
1933  fHistoDoubleCountTrueClusterGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt(),(Double_t)0,fWeightJetJetMC);
1934  FillMultipleCountMap(fMapMultipleCountTrueClusterGammas,motherLab);
1935  }
1936  }
1937  Int_t grandMotherLab = ((AliAODMCParticle*) AODMCTrackArray->At(motherLab))->GetMother();
1938  if (grandMotherLab > -1){
1939  if (TruePhotonCandidate->IsLargestComponentElectron() && TruePhotonCandidate->IsConversion()){
1940  if (CheckVectorForDoubleCount(fVectorDoubleCountTrueClusterGammas,grandMotherLab)){
1941  fHistoDoubleCountTrueClusterGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt(),(Double_t)1,fWeightJetJetMC);
1942  FillMultipleCountMap(fMapMultipleCountTrueClusterGammas,grandMotherLab);
1943  }
1944  }
1945  }
1946  }
1947 
1948  return;
1949 }
1950 
1951 //________________________________________________________________________
1952 void AliAnalysisTaskHeavyNeutralMesonToGG::ProcessPhotonCandidates(){
1953 
1954  Int_t nV0 = 0;
1955  TList *GammaCandidatesStepOne = new TList();
1956  TList *GammaCandidatesStepTwo = new TList();
1957  // Loop over Photon Candidates allocated by ReaderV1
1958  for(Int_t i = 0; i < fReaderGammas->GetEntriesFast(); i++){
1959  AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) fReaderGammas->At(i);
1960  if(!PhotonCandidate) continue;
1961  fIsFromDesiredHeader = kTRUE;
1962  if(fIsMC>0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0){
1963  Int_t isPosFromMBHeader = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCEvent, fInputEvent);
1964  if(isPosFromMBHeader == 0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1965  Int_t isNegFromMBHeader = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCEvent, fInputEvent);
1966  if(isNegFromMBHeader == 0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1967  if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromDesiredHeader = kFALSE;
1968  }
1969 
1970  if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelected(PhotonCandidate,fInputEvent)) continue;
1971  if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(PhotonCandidate->GetPhotonPhi(),fEventPlaneAngle)) continue;
1972  if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseElecSharingCut() &&
1973  !((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){
1974  fGammaCandidates->Add(PhotonCandidate); // if no second loop is required add to events good gammas
1975 
1976  if(fIsFromDesiredHeader){
1977  if(!fDoLightOutput) fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt(),fWeightJetJetMC);
1978  }
1979  if(fIsMC>0){
1980  if(fInputEvent->IsA()==AliESDEvent::Class())
1981  ProcessTruePhotonCandidates(PhotonCandidate);
1982  if(fInputEvent->IsA()==AliAODEvent::Class())
1983  ProcessTruePhotonCandidatesAOD(PhotonCandidate);
1984  }
1985  }else if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseElecSharingCut()){ // if Shared Electron cut is enabled, Fill array, add to step one
1986  ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->FillElectonLabelArray(PhotonCandidate,nV0);
1987  nV0++;
1988  GammaCandidatesStepOne->Add(PhotonCandidate);
1989  }else if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseElecSharingCut() &&
1990  ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){ // shared electron is disabled, step one not needed -> step two
1991  GammaCandidatesStepTwo->Add(PhotonCandidate);
1992  }
1993  }
1994 
1995  if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseElecSharingCut()){
1996  for(Int_t i = 0;i<GammaCandidatesStepOne->GetEntries();i++){
1997  AliAODConversionPhoton *PhotonCandidate= (AliAODConversionPhoton*) GammaCandidatesStepOne->At(i);
1998  if(!PhotonCandidate) continue;
1999  fIsFromDesiredHeader = kTRUE;
2000  if(fMCEvent && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0){
2001  Int_t isPosFromMBHeader = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCEvent, fInputEvent);
2002  Int_t isNegFromMBHeader = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCEvent, fInputEvent);
2003  if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromDesiredHeader = kFALSE;
2004  }
2005  if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->RejectSharedElectronV0s(PhotonCandidate,i,GammaCandidatesStepOne->GetEntries())) continue;
2006  if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){ // To Colse v0s cut diabled, step two not needed
2007  fGammaCandidates->Add(PhotonCandidate);
2008  if(fIsFromDesiredHeader){
2009  if(!fDoLightOutput) fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt(),fWeightJetJetMC);
2010  }
2011  if(fIsMC>0){
2012  if(fInputEvent->IsA()==AliESDEvent::Class())
2013  ProcessTruePhotonCandidates(PhotonCandidate);
2014  if(fInputEvent->IsA()==AliAODEvent::Class())
2015  ProcessTruePhotonCandidatesAOD(PhotonCandidate);
2016  }
2017  } else GammaCandidatesStepTwo->Add(PhotonCandidate); // Close v0s cut enabled -> add to list two
2018  }
2019  }
2020 
2021  if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){
2022  for(Int_t i = 0;i<GammaCandidatesStepTwo->GetEntries();i++){
2023  AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) GammaCandidatesStepTwo->At(i);
2024  if(!PhotonCandidate) continue;
2025  fIsFromDesiredHeader = kTRUE;
2026  if(fMCEvent && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0){
2027  Int_t isPosFromMBHeader = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCEvent, fInputEvent);
2028  Int_t isNegFromMBHeader = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCEvent, fInputEvent);
2029  if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromDesiredHeader = kFALSE;
2030  }
2031  if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->RejectToCloseV0s(PhotonCandidate,GammaCandidatesStepTwo,i)) continue;
2032  fGammaCandidates->Add(PhotonCandidate); // Add gamma to current cut TList
2033  if(fIsFromDesiredHeader){
2034  if(!fDoLightOutput) fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt(),fWeightJetJetMC);
2035  }
2036  if(fIsMC>0){
2037  if(fInputEvent->IsA()==AliESDEvent::Class())
2038  ProcessTruePhotonCandidates(PhotonCandidate);
2039  if(fInputEvent->IsA()==AliAODEvent::Class())
2040  ProcessTruePhotonCandidatesAOD(PhotonCandidate);
2041  }
2042  }
2043  }
2044 
2045  delete GammaCandidatesStepOne;
2046  GammaCandidatesStepOne = 0x0;
2047  delete GammaCandidatesStepTwo;
2048  GammaCandidatesStepTwo = 0x0;
2049 
2050 }
2051 
2052 //________________________________________________________________________
2053 void AliAnalysisTaskHeavyNeutralMesonToGG::ProcessTruePhotonCandidatesAOD(AliAODConversionPhoton *TruePhotonCandidate){
2054 
2055  const AliVVertex* primVtxMC = fMCEvent->GetPrimaryVertex();
2056  Double_t mcProdVtxX = primVtxMC->GetX();
2057  Double_t mcProdVtxY = primVtxMC->GetY();
2058  Double_t mcProdVtxZ = primVtxMC->GetZ();
2059 
2060  TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
2061  if (AODMCTrackArray == NULL) return;
2062  AliAODMCParticle *posDaughter = (AliAODMCParticle*) AODMCTrackArray->At(TruePhotonCandidate->GetMCLabelPositive());
2063  AliAODMCParticle *negDaughter = (AliAODMCParticle*) AODMCTrackArray->At(TruePhotonCandidate->GetMCLabelNegative());
2064 
2065  if(posDaughter == NULL || negDaughter == NULL) return; // One particle does not exist
2066  Int_t pdgCode[2] = {TMath::Abs(posDaughter->GetPdgCode()),TMath::Abs(negDaughter->GetPdgCode())};
2067 
2068  if(posDaughter->GetMother() != negDaughter->GetMother()){
2069  return;
2070  } else if(posDaughter->GetMother() == -1){
2071  return;
2072  }
2073 
2074  if(pdgCode[0]!=11 || pdgCode[1]!=11){
2075  return; //One Particle is not a electron
2076  }
2077 
2078  if(posDaughter->GetPdgCode()==negDaughter->GetPdgCode()){
2079  return; // Same Charge
2080  }
2081 
2082  AliAODMCParticle *Photon = (AliAODMCParticle*) AODMCTrackArray->At(posDaughter->GetMother());
2083  if(Photon->GetPdgCode() != 22){
2084  return; // Mother is no Photon
2085  }
2086 
2087  if(((posDaughter->GetMCProcessCode())) != 5 || ((negDaughter->GetMCProcessCode())) != 5){
2088  return;// check if the daughters come from a conversion
2089  }
2090  // STILL A BUG IN ALIROOT >>8 HAS TPO BE REMOVED AFTER FIX
2091 
2092  // True Photon
2093  if(fIsFromDesiredHeader){
2094  if(!fDoLightOutput) fHistoTrueConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt(),fWeightJetJetMC);
2095  if (CheckVectorForDoubleCount(fVectorDoubleCountTrueConvGammas,posDaughter->GetMother())){
2096  if(!fDoLightOutput) fHistoDoubleCountTrueConvGammaRPt[fiCut]->Fill(TruePhotonCandidate->GetConversionRadius(),TruePhotonCandidate->Pt(),fWeightJetJetMC);
2097  FillMultipleCountMap(fMapMultipleCountTrueConvGammas,posDaughter->GetMother());
2098  }
2099  }
2100 
2101  Bool_t isPrimary = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsConversionPrimaryAOD(fInputEvent, Photon, mcProdVtxX, mcProdVtxY, mcProdVtxZ);
2102  if(isPrimary){
2103  // Count just primary MC Gammas as true --> For Ratio esdtruegamma / mcconvgamma
2104  if(fIsFromDesiredHeader){
2105  if(!fDoLightOutput){
2106  fHistoTruePrimaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt(),fWeightJetJetMC);
2107  }
2108  }
2109  // (Not Filled for i6, Extra Signal Gamma (parambox) are secondary)
2110  }
2111  TruePhotonCandidate->SetIsTrueConvertedPhoton();
2112  return;
2113 }
2114 
2115 //________________________________________________________________________
2116 void AliAnalysisTaskHeavyNeutralMesonToGG::ProcessTruePhotonCandidates(AliAODConversionPhoton *TruePhotonCandidate){
2117 
2118  const AliVVertex* primVtxMC = fMCEvent->GetPrimaryVertex();
2119  Double_t mcProdVtxX = primVtxMC->GetX();
2120  Double_t mcProdVtxY = primVtxMC->GetY();
2121  Double_t mcProdVtxZ = primVtxMC->GetZ();
2122 
2123  // Process True Photons
2124  TParticle *posDaughter = TruePhotonCandidate->GetPositiveMCDaughter(fMCEvent);
2125  TParticle *negDaughter = TruePhotonCandidate->GetNegativeMCDaughter(fMCEvent);
2126 
2127  if(posDaughter == NULL || negDaughter == NULL) return; // One particle does not exist
2128  Int_t pdgCode[2] = {TMath::Abs(posDaughter->GetPdgCode()),TMath::Abs(negDaughter->GetPdgCode())};
2129  if(posDaughter->GetMother(0) != negDaughter->GetMother(0)){
2130  return;
2131  }
2132  else if(posDaughter->GetMother(0) == -1){
2133  return;
2134  }
2135 
2136  if(pdgCode[0]!=11 || pdgCode[1]!=11) return; //One Particle is not a electron
2137 
2138  if(posDaughter->GetPdgCode()==negDaughter->GetPdgCode()) return; // Same Charge
2139 
2140  TParticle *Photon = TruePhotonCandidate->GetMCParticle(fMCEvent);
2141 
2142  if(Photon->GetPdgCode() != 22){
2143  return; // Mother is no Photon
2144  }
2145 
2146  if(posDaughter->GetUniqueID() != 5 || negDaughter->GetUniqueID() !=5) return;// check if the daughters come from a conversion
2147 
2148  // True Photon
2149  if(fIsFromDesiredHeader){
2150  if(!fDoLightOutput) fHistoTrueConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt(),fWeightJetJetMC);
2151  if (CheckVectorForDoubleCount(fVectorDoubleCountTrueConvGammas,posDaughter->GetMother(0))){
2152  if(!fDoLightOutput) fHistoDoubleCountTrueConvGammaRPt[fiCut]->Fill(TruePhotonCandidate->GetConversionRadius(),TruePhotonCandidate->Pt(),fWeightJetJetMC);
2153  FillMultipleCountMap(fMapMultipleCountTrueConvGammas,posDaughter->GetMother(0));
2154  }
2155  }
2156  Bool_t isPrimary = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsConversionPrimaryESD( fMCEvent, posDaughter->GetMother(0), mcProdVtxX, mcProdVtxY, mcProdVtxZ);
2157  if(isPrimary){
2158  // Count just primary MC Gammas as true --> For Ratio esdtruegamma / mcconvgamma
2159  if(fIsFromDesiredHeader){
2160  if(!fDoLightOutput){
2161  fHistoTruePrimaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt(),fWeightJetJetMC);
2162  }
2163  }
2164  // (Not Filled for i6, Extra Signal Gamma (parambox) are secondary)
2165  }
2166  TruePhotonCandidate->SetIsTrueConvertedPhoton();
2167  return;
2168 }
2169 //________________________________________________________________________
2170 void AliAnalysisTaskHeavyNeutralMesonToGG::ProcessAODMCParticles(){
2171 
2172  const AliVVertex* primVtxMC = fMCEvent->GetPrimaryVertex();
2173  Double_t mcProdVtxX = primVtxMC->GetX();
2174  Double_t mcProdVtxY = primVtxMC->GetY();
2175  Double_t mcProdVtxZ = primVtxMC->GetZ();
2176 
2177  TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
2178  if (AODMCTrackArray == NULL) return;
2179 
2180  // Loop over all primary MC particle
2181  for(Long_t i = 0; i < AODMCTrackArray->GetEntriesFast(); i++) {
2182 
2183  AliAODMCParticle* particle = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(i));
2184  if (!particle) continue;
2185 
2186  Bool_t isPrimary = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsConversionPrimaryAOD(fInputEvent, particle, mcProdVtxX, mcProdVtxY, mcProdVtxZ);
2187  if (isPrimary) {
2188 
2189  Int_t isMCFromMBHeader = -1;
2190  if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0){
2191  isMCFromMBHeader = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCEvent, fInputEvent);
2192  if(isMCFromMBHeader == 0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
2193  }
2194 
2195  if (fMesonRecoMode < 2){
2196  if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(particle->Phi(),fEventPlaneAngle,kFALSE)) continue;
2197  }
2198 
2199  Double_t mesonY = 1.e30;
2200  Double_t ratio = 0;
2201  if (particle->E() != TMath::Abs(particle->Pz())){
2202  ratio = (particle->E()+particle->Pz()) / (particle->E()-particle->Pz());
2203  }
2204  if( !(ratio <= 0) ){
2205  mesonY = particle->Y()-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift();
2206  }
2207 
2208  // check neutral mesons
2209  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelectedAODMC(particle,AODMCTrackArray,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift())){
2210  AliAODMCParticle* daughter0 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(particle->GetDaughter(0)));
2211  AliAODMCParticle* daughter1 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(particle->GetDaughter(1)));
2212  Float_t weighted= 1;
2213  if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCEvent, fInputEvent)){
2214  if (particle->Pt()>0.005){
2215  weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(i, 0x0, fInputEvent);
2216  }
2217  }
2218  Double_t alpha = -10;
2219 
2220  if(particle->GetPdgCode() == fMesonPDG){
2221  alpha = (daughter0->E() - daughter1->E())/(daughter0->E() + daughter1->E());
2222  if (fHistoMCMesonPt[fiCut])fHistoMCMesonPt[fiCut]->Fill(particle->Pt(),weighted*fWeightJetJetMC); // All MC Meson
2223  if (fHistoMCMesonWOWeightPt[fiCut]) fHistoMCMesonWOWeightPt[fiCut]->Fill(particle->Pt(),fWeightJetJetMC);
2224  if (fIsMC > 1 && fHistoMCMesonWOEvtWeightPt[fiCut])fHistoMCMesonWOEvtWeightPt[fiCut]->Fill(particle->Pt());
2225  if (fDoMesonQA > 0){
2226  if (fHistoMCMesonPtY[fiCut])fHistoMCMesonPtY[fiCut]->Fill(particle->Pt(),mesonY,weighted*fWeightJetJetMC);
2227  if (fHistoMCMesonPtAlpha[fiCut])fHistoMCMesonPtAlpha[fiCut]->Fill(particle->Pt(),alpha,fWeightJetJetMC);
2228  if (fIsMC == 2 && fHistoMCMesonPtJetPt[fiCut])fHistoMCMesonPtJetPt[fiCut]->Fill(particle->Pt(),((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetMaxPtJet(),fWeightJetJetMC);
2229  }
2230 
2231  // Check the acceptance for both gammas
2232  if (fMesonRecoMode == 0){
2233  if( ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(daughter0,AODMCTrackArray,kFALSE) &&
2234  ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(daughter1,AODMCTrackArray,kFALSE) &&
2235  ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter0->Phi(),fEventPlaneAngle,kFALSE) &&
2236  ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter1->Phi(),fEventPlaneAngle,kFALSE)){
2237  // check acceptance of clusters as well, true if one of them points into the Calo acceptance
2238  if (fHistoMCMesonInAccPt[fiCut])fHistoMCMesonInAccPt[fiCut]->Fill(particle->Pt(),weighted*fWeightJetJetMC); // MC Meson with gamma in acc
2239  if (fHistoMCMesonWOWeightInAccPt[fiCut])fHistoMCMesonWOWeightInAccPt[fiCut]->Fill(particle->Pt(),fWeightJetJetMC); // MC Meson with gamma in acc wo weight
2240  if(fIsMC > 1 && fHistoMCMesonWOEvtWeightInAccPt[fiCut])fHistoMCMesonWOEvtWeightInAccPt[fiCut]->Fill(particle->Pt()); // MC Meson with gamma in acc wo any weight
2241  }
2242  } else if (fMesonRecoMode == 1){
2243  if( ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(daughter0,AODMCTrackArray,kFALSE) &&
2244  ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(daughter1,AODMCTrackArray,kFALSE) &&
2245  ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter0->Phi(),fEventPlaneAngle,kFALSE) &&
2246  ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter1->Phi(),fEventPlaneAngle,kFALSE)){
2247  // check acceptance of clusters as well, true if one of them points into the Calo acceptance
2248  if( ((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->ClusterIsSelectedAODMC(daughter0,AODMCTrackArray) ||
2249  ((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->ClusterIsSelectedAODMC(daughter1,AODMCTrackArray) ){
2250  if (fHistoMCMesonInAccPt[fiCut])fHistoMCMesonInAccPt[fiCut]->Fill(particle->Pt(),weighted*fWeightJetJetMC); // MC Meson with gamma in acc
2251  if (fHistoMCMesonWOWeightInAccPt[fiCut])fHistoMCMesonWOWeightInAccPt[fiCut]->Fill(particle->Pt(),fWeightJetJetMC); // MC Meson with gamma in acc wo weight
2252  if(fIsMC > 1 && fHistoMCMesonWOEvtWeightInAccPt[fiCut])fHistoMCMesonWOEvtWeightInAccPt[fiCut]->Fill(particle->Pt()); // MC Meson with gamma in acc wo any weight
2253  }
2254  }
2255  } else if (fMesonRecoMode == 2){
2256  if (((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->ClusterIsSelectedAODMC(daughter0,AODMCTrackArray) &&
2257  ((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->ClusterIsSelectedAODMC(daughter1,AODMCTrackArray) ){
2258  if (fHistoMCMesonInAccPt[fiCut])fHistoMCMesonInAccPt[fiCut]->Fill(particle->Pt(),weighted*fWeightJetJetMC); // MC Meson with gamma in acc
2259  if (fHistoMCMesonWOWeightInAccPt[fiCut])fHistoMCMesonWOWeightInAccPt[fiCut]->Fill(particle->Pt(),fWeightJetJetMC); // MC Meson with gamma in acc wo weight
2260  if(fIsMC > 1 && fHistoMCMesonWOEvtWeightInAccPt[fiCut])fHistoMCMesonWOEvtWeightInAccPt[fiCut]->Fill(particle->Pt()); // MC Meson with gamma in acc wo any weight
2261  }
2262  }
2263  }
2264  }
2265  }
2266  }
2267 }
2268 
2269 //________________________________________________________________________
2270 void AliAnalysisTaskHeavyNeutralMesonToGG::ProcessMCParticles(){
2271 
2272  const AliVVertex* primVtxMC = fMCEvent->GetPrimaryVertex();
2273  Double_t mcProdVtxX = primVtxMC->GetX();
2274  Double_t mcProdVtxY = primVtxMC->GetY();
2275  Double_t mcProdVtxZ = primVtxMC->GetZ();
2276  // cout << mcProdVtxX <<"\t" << mcProdVtxY << "\t" << mcProdVtxZ << endl;
2277 
2278  // Loop over all primary MC particle
2279  for(Long_t i = 0; i < fMCEvent->GetNumberOfTracks(); i++) {
2280  if (((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsConversionPrimaryESD( fMCEvent, i, mcProdVtxX, mcProdVtxY, mcProdVtxZ)){
2281  // fill primary histograms
2282  TParticle* particle = (TParticle *)fMCEvent->Particle(i);
2283  if (!particle) continue;
2284 
2285  Int_t isMCFromMBHeader = -1;
2286  if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0){
2287  isMCFromMBHeader = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCEvent, fInputEvent);
2288  if(isMCFromMBHeader == 0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
2289  }
2290 
2291  if (fMesonRecoMode < 2){
2292  if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(particle->Phi(),fEventPlaneAngle,kFALSE)) continue;
2293  }
2294 
2295  // Fill histograms for other particles
2296  Double_t mesonY = 1.e30;
2297  Double_t ratio = 0;
2298  if (particle->Energy() != TMath::Abs(particle->Pz())){
2299  ratio = (particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz());
2300  }
2301  if( !(ratio <= 0) ){
2302  mesonY = particle->Y()-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift();
2303  }
2304 
2305  // check neutral mesons
2306  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelectedMC(particle,fMCEvent,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift())){
2307  TParticle* daughter0 = (TParticle*)fMCEvent->Particle(particle->GetFirstDaughter());
2308  TParticle* daughter1 = (TParticle*)fMCEvent->Particle(particle->GetLastDaughter());
2309 
2310  Float_t weighted= 1;
2311  if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCEvent, fInputEvent)){
2312  if (particle->Pt()>0.005){
2313  weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(i, fMCEvent, fInputEvent);
2314  }
2315  }
2316 
2317  Double_t alpha = -10;
2318  if(particle->GetPdgCode() == fMesonPDG){
2319  alpha = (daughter0->Energy() - daughter1->Energy())/(daughter0->Energy() + daughter1->Energy());
2320  if (fHistoMCMesonPt[fiCut])fHistoMCMesonPt[fiCut]->Fill(particle->Pt(),weighted*fWeightJetJetMC); // All MC Meson
2321  if (fHistoMCMesonWOWeightPt[fiCut]) fHistoMCMesonWOWeightPt[fiCut]->Fill(particle->Pt(),fWeightJetJetMC);
2322  if (fIsMC > 1 && fHistoMCMesonWOEvtWeightPt[fiCut])fHistoMCMesonWOEvtWeightPt[fiCut]->Fill(particle->Pt());
2323  if (fDoMesonQA > 0){
2324  if (fHistoMCMesonPtY[fiCut])fHistoMCMesonPtY[fiCut]->Fill(particle->Pt(),mesonY,weighted*fWeightJetJetMC); // All MC Meson
2325  if (fHistoMCMesonPtAlpha[fiCut])fHistoMCMesonPtAlpha[fiCut]->Fill(particle->Pt(),alpha,fWeightJetJetMC); // All MC Meson
2326  if (fIsMC == 2 && fHistoMCMesonPtJetPt[fiCut])fHistoMCMesonPtJetPt[fiCut]->Fill(particle->Pt(),((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetMaxPtJet(),fWeightJetJetMC);
2327  }
2328 
2329  // Check the acceptance for both gammas & whether they are counted as primaries as well
2330  Bool_t kDaughter0IsPrim = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsConversionPrimaryESD( fMCEvent, particle->GetFirstDaughter(), mcProdVtxX, mcProdVtxY, mcProdVtxZ);
2331  Bool_t kDaughter1IsPrim = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsConversionPrimaryESD( fMCEvent, particle->GetLastDaughter(), mcProdVtxX, mcProdVtxY, mcProdVtxZ);
2332  if( kDaughter0IsPrim && kDaughter1IsPrim ){
2333  if (fMesonRecoMode == 0){
2334  if( ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(daughter0,fMCEvent,kFALSE) &&
2335  ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(daughter1,fMCEvent,kFALSE) &&
2336  ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter0->Phi(),fEventPlaneAngle,kFALSE) &&
2337  ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter1->Phi(),fEventPlaneAngle,kFALSE)){
2338  // check acceptance of clusters as well, true if one of them points into the Calo acceptance
2339  if (fHistoMCMesonInAccPt[fiCut])fHistoMCMesonInAccPt[fiCut]->Fill(particle->Pt(),weighted*fWeightJetJetMC); // MC Meson with gamma in acc
2340  if (fHistoMCMesonWOWeightInAccPt[fiCut])fHistoMCMesonWOWeightInAccPt[fiCut]->Fill(particle->Pt(),fWeightJetJetMC); // MC Meson with gamma in acc wo weighting
2341  if(fIsMC > 1 && fHistoMCMesonWOEvtWeightInAccPt[fiCut])fHistoMCMesonWOEvtWeightInAccPt[fiCut]->Fill(particle->Pt()); // MC Meson with gamma in acc wo any weight
2342  }
2343  } else if (fMesonRecoMode == 1){
2344  if (((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(daughter0,fMCEvent,kFALSE) &&
2345  ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(daughter1,fMCEvent,kFALSE) &&
2346  ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter0->Phi(),fEventPlaneAngle,kFALSE) &&
2347  ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter1->Phi(),fEventPlaneAngle,kFALSE)){
2348  // check acceptance of clusters as well, true if one of them points into the Calo acceptance
2349  if (((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->ClusterIsSelectedMC(daughter0,fMCEvent) ||
2350  ((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->ClusterIsSelectedMC(daughter1,fMCEvent) ){
2351  if (fHistoMCMesonInAccPt[fiCut])fHistoMCMesonInAccPt[fiCut]->Fill(particle->Pt(),weighted*fWeightJetJetMC); // MC Meson with gamma in acc
2352  if (fHistoMCMesonWOWeightInAccPt[fiCut])fHistoMCMesonWOWeightInAccPt[fiCut]->Fill(particle->Pt(),fWeightJetJetMC); // MC Meson with gamma in acc wo weighting
2353  if(fIsMC > 1 && fHistoMCMesonWOEvtWeightInAccPt[fiCut])fHistoMCMesonWOEvtWeightInAccPt[fiCut]->Fill(particle->Pt()); // MC Meson with gamma in acc wo any weight
2354  }
2355  }
2356  } else if (fMesonRecoMode == 2){
2357  if (((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->ClusterIsSelectedMC(daughter0,fMCEvent) &&
2358  ((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->ClusterIsSelectedMC(daughter1,fMCEvent) ){
2359  if (fHistoMCMesonInAccPt[fiCut])fHistoMCMesonInAccPt[fiCut]->Fill(particle->Pt(),weighted*fWeightJetJetMC); // MC Meson with gamma in acc
2360  if (fHistoMCMesonWOWeightInAccPt[fiCut])fHistoMCMesonWOWeightInAccPt[fiCut]->Fill(particle->Pt(),fWeightJetJetMC); // MC Meson with gamma in acc wo weighting
2361  if(fIsMC > 1 && fHistoMCMesonWOEvtWeightInAccPt[fiCut])fHistoMCMesonWOEvtWeightInAccPt[fiCut]->Fill(particle->Pt()); // MC Meson with gamma in acc wo any weight
2362  }
2363  }
2364  }
2365  }
2366  }
2367  }
2368  }
2369 }
2370 
2371 //________________________________________________________________________
2372 void AliAnalysisTaskHeavyNeutralMesonToGG::CalculateMesonCandidates(){
2373  TClonesArray * arrClustersMesonCand = NULL;
2374  if(fCorrTaskSetting.CompareTo("") && fMesonRecoMode > 0)
2375  arrClustersMesonCand = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(Form("%sClustersBranch",fCorrTaskSetting.Data())));
2376 
2377  if (fMesonRecoMode == 0){
2378  if(fGammaCandidates->GetEntries()>1){
2379  for(Int_t firstGammaIndex=0;firstGammaIndex<fGammaCandidates->GetEntries()-1;firstGammaIndex++){
2380  AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(firstGammaIndex));
2381  if (gamma0==NULL) continue;
2382  for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fGammaCandidates->GetEntries();secondGammaIndex++){
2383  AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(secondGammaIndex));
2384  //Check for same Electron ID
2385  if (gamma1==NULL) continue;
2386  if(gamma0->GetTrackLabelPositive() == gamma1->GetTrackLabelPositive() ||
2387  gamma0->GetTrackLabelNegative() == gamma1->GetTrackLabelNegative() ||
2388  gamma0->GetTrackLabelNegative() == gamma1->GetTrackLabelPositive() ||
2389  gamma0->GetTrackLabelPositive() == gamma1->GetTrackLabelNegative() ) continue;
2390 
2391  AliAODConversionMother *mesonCand = new AliAODConversionMother(gamma0,gamma1);
2392  mesonCand->SetLabels(firstGammaIndex,secondGammaIndex);
2393  mesonCand->CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
2394 
2395  if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(mesonCand,kTRUE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift()))){
2396  if (fHistoMotherInvMassPt[fiCut]) fHistoMotherInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
2397 
2398  if (fDoMesonQA > 0){
2399  if (mesonCand->M() > fMesonInvMassWindow[0] && mesonCand->M() < fMesonInvMassWindow[1]){
2400  if (fHistoMotherMesonPtY[fiCut]) fHistoMotherMesonPtY[fiCut]->Fill(mesonCand->Pt(),mesonCand->Rapidity()-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift(),fWeightJetJetMC);
2401  if (fHistoMotherMesonPtAlpha[fiCut]) fHistoMotherMesonPtAlpha[fiCut]->Fill(mesonCand->Pt(),mesonCand->GetAlpha(),fWeightJetJetMC);
2402  if (fHistoMotherMesonPtOpenAngle[fiCut]) fHistoMotherMesonPtOpenAngle[fiCut]->Fill(mesonCand->Pt(),mesonCand->GetOpeningAngle(),fWeightJetJetMC);
2403  }
2404  }
2405  if(fDoTHnSparse && ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->DoBGCalculation()){
2406  Int_t zbin = 0;
2407  Int_t mbin = 0;
2408  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->BackgroundHandlerType() == 0){
2409  zbin = fBGHandler[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
2410  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2411  mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
2412  }else {
2413  mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries());
2414  }
2415  }else{
2416  zbin = fBGHandlerRP[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
2417  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2418  mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
2419  }else {
2420  mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries());
2421  }
2422  }
2423 
2424  Double_t sparesFill[4] = {mesonCand->M(),mesonCand->Pt(),(Double_t)zbin,(Double_t)mbin};
2425  if (fSparseMotherInvMassPtZM[fiCut]) fSparseMotherInvMassPtZM[fiCut]->Fill(sparesFill,1);
2426  }
2427 
2428 
2429  if( fIsMC > 0 ){
2430  if(fInputEvent->IsA()==AliESDEvent::Class())
2431  ProcessTrueMesonCandidatesConv(mesonCand, gamma0, gamma1);
2432  if(fInputEvent->IsA()==AliAODEvent::Class())
2433  ProcessTrueMesonCandidatesConvAOD(mesonCand, gamma0, gamma1);
2434  }
2435  }
2436  delete mesonCand;
2437  mesonCand=0x0;
2438  }
2439  }
2440  }
2441  } else if (fMesonRecoMode == 1){
2442  if(fGammaCandidates->GetEntries()>0){
2443  for(Int_t firstGammaIndex=0;firstGammaIndex<fGammaCandidates->GetEntries();firstGammaIndex++){
2444  AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(firstGammaIndex));
2445  if (gamma0==NULL) continue;
2446 
2447  for(Int_t secondGammaIndex=0;secondGammaIndex<fClusterCandidates->GetEntries();secondGammaIndex++){
2448  Bool_t matched = kFALSE;
2449  AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fClusterCandidates->At(secondGammaIndex));
2450  if (gamma1==NULL) continue;
2451 
2452  if (gamma1->GetIsCaloPhoton()){
2453  AliVCluster* cluster = NULL;
2454  if(fInputEvent->IsA()==AliESDEvent::Class()){
2455  if(arrClustersMesonCand)
2456  cluster = new AliESDCaloCluster(*(AliESDCaloCluster*)arrClustersMesonCand->At(gamma1->GetCaloClusterRef()));
2457  else
2458  cluster = fInputEvent->GetCaloCluster(gamma1->GetCaloClusterRef());
2459  } else if(fInputEvent->IsA()==AliAODEvent::Class()){
2460  if(arrClustersMesonCand)
2461  cluster = new AliAODCaloCluster(*(AliAODCaloCluster*)arrClustersMesonCand->At(gamma1->GetCaloClusterRef()));
2462  else
2463  cluster = fInputEvent->GetCaloCluster(gamma1->GetCaloClusterRef());
2464  }
2465 
2466  matched = ((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->MatchConvPhotonToCluster(gamma0,cluster, fInputEvent, fWeightJetJetMC);
2467  if(arrClustersMesonCand) delete cluster;
2468  }
2469 
2470  AliAODConversionMother *mesonCand = new AliAODConversionMother(gamma0,gamma1);
2471  mesonCand->SetLabels(firstGammaIndex,secondGammaIndex);
2472 
2473  if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(mesonCand,kTRUE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift()))){
2474  if (matched){
2475  if(fHistoMotherMatchedInvMassPt[fiCut]) fHistoMotherMatchedInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),fWeightJetJetMC);
2476  }else {
2477  if (fHistoMotherInvMassPt[fiCut]) fHistoMotherInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),fWeightJetJetMC);
2478  }
2479 
2480  // fill new histograms
2481  if (!matched){
2482  if (fDoMesonQA > 0){
2483  if (mesonCand->M() > fMesonInvMassWindow[0] && mesonCand->M() < fMesonInvMassWindow[1]){
2484  if (fHistoMotherMesonPtY[fiCut]) fHistoMotherMesonPtY[fiCut]->Fill(mesonCand->Pt(),mesonCand->Rapidity()-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift(),fWeightJetJetMC);
2485  if (fHistoMotherMesonPtAlpha[fiCut]) fHistoMotherMesonPtAlpha[fiCut]->Fill(mesonCand->Pt(),mesonCand->GetAlpha(),fWeightJetJetMC);
2486  if (fHistoMotherMesonPtOpenAngle[fiCut]) fHistoMotherMesonPtOpenAngle[fiCut]->Fill(mesonCand->Pt(),mesonCand->GetOpeningAngle(),fWeightJetJetMC);
2487  if (fHistoMotherMesonConvPhotonEtaPhi[fiCut]) fHistoMotherMesonConvPhotonEtaPhi[fiCut]->Fill(gamma0->GetPhotonPhi(), gamma0->GetPhotonEta(),fWeightJetJetMC);
2488  }
2489  }
2490  if(fDoTHnSparse && ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->DoBGCalculation()){
2491  Int_t zbin = 0;
2492  Int_t mbin = 0;
2493 
2494  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->BackgroundHandlerType() == 0){
2495  zbin = fBGHandler[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
2496  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2497  mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
2498  }else {
2499  mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries());
2500  }
2501  }else{
2502  zbin = fBGHandlerRP[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
2503  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2504  mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
2505  }else {
2506  mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries());
2507  }
2508  }
2509  Double_t sparesFill[4] = {mesonCand->M(),mesonCand->Pt(),(Double_t)zbin,(Double_t)mbin};
2510  if (fSparseMotherInvMassPtZM[fiCut]) fSparseMotherInvMassPtZM[fiCut]->Fill(sparesFill,1);
2511  }
2512  }
2513 
2514  if(fIsMC>0){
2515  if(fInputEvent->IsA()==AliESDEvent::Class())
2516  ProcessTrueMesonCandidatesConvCalo(mesonCand, gamma0, gamma1, matched);
2517  if(fInputEvent->IsA()==AliAODEvent::Class())
2518  ProcessTrueMesonCandidatesConvCaloAOD(mesonCand, gamma0, gamma1, matched);
2519  }
2520  }
2521  delete mesonCand;
2522  mesonCand=0x0;
2523  }
2524  }
2525  }
2526  } else if (fMesonRecoMode == 2){
2527  if(fClusterCandidates->GetEntries()>0){
2528 
2529  for(Int_t firstGammaIndex=0;firstGammaIndex<fClusterCandidates->GetEntries();firstGammaIndex++){
2530  AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fClusterCandidates->At(firstGammaIndex));
2531  if (gamma0==NULL) continue;
2532  for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fClusterCandidates->GetEntries();secondGammaIndex++){
2533  AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fClusterCandidates->At(secondGammaIndex));
2534  if (gamma1==NULL) continue;
2535 
2536  AliAODConversionMother *mesonCand = new AliAODConversionMother(gamma0,gamma1);
2537  mesonCand->SetLabels(firstGammaIndex,secondGammaIndex);
2538 
2539  if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(mesonCand, kTRUE, ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift(), gamma0->GetLeadingCellID(), gamma1->GetLeadingCellID()))){
2540  if (fHistoMotherInvMassPt[fiCut]) fHistoMotherInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
2541  // fill new histograms
2542  if (fDoMesonQA > 0 && fDoMesonQA < 3){
2543  if ( mesonCand->M() > fMesonInvMassWindow[0] && mesonCand->M() < fMesonInvMassWindow[1]){
2544  fHistoMotherMesonPtY[fiCut]->Fill(mesonCand->Pt(),mesonCand->Rapidity()-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift(), fWeightJetJetMC);
2545  fHistoMotherMesonPtAlpha[fiCut]->Fill(mesonCand->Pt(),TMath::Abs(mesonCand->GetAlpha()), fWeightJetJetMC);
2546  fHistoMotherMesonPtOpenAngle[fiCut]->Fill(mesonCand->Pt(),mesonCand->GetOpeningAngle(), fWeightJetJetMC);
2547  }
2548  }
2549  if(fDoTHnSparse && ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->DoBGCalculation()){
2550  Int_t zbin = 0;
2551  Int_t mbin = 0;
2552 
2553  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->BackgroundHandlerType() == 0){
2554  zbin = fBGHandler[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
2555  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2556  mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
2557  } else {
2558  mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fClusterCandidates->GetEntries());
2559  }
2560  }
2561  Double_t sparesFill[4] = {mesonCand->M(),mesonCand->Pt(),(Double_t)zbin,(Double_t)mbin};
2562  if (fSparseMotherInvMassPtZM[fiCut]) fSparseMotherInvMassPtZM[fiCut]->Fill(sparesFill,1);
2563  }
2564 
2565  if(fIsMC> 0){
2566  if(fInputEvent->IsA()==AliESDEvent::Class())
2567  ProcessTrueMesonCandidatesCalo(mesonCand,gamma0,gamma1);
2568  if(fInputEvent->IsA()==AliAODEvent::Class())
2569  ProcessTrueMesonCandidatesCaloAOD(mesonCand,gamma0,gamma1);
2570  }
2571  }
2572  delete mesonCand;
2573  mesonCand=0x0;
2574  }
2575  }
2576  }
2577  }
2578 }
2579 
2580 //______________________________________________________________________
2581 void AliAnalysisTaskHeavyNeutralMesonToGG::ProcessTrueMesonCandidatesConvCalo( AliAODConversionMother *mesonCand,
2582  AliAODConversionPhoton *TrueGammaCandidate0,
2583  AliAODConversionPhoton *TrueGammaCandidate1,
2584  Bool_t matched )
2585 {
2586  // obtain MC vertex
2587  const AliVVertex* primVtxMC = fMCEvent->GetPrimaryVertex();
2588  Double_t mcProdVtxX = primVtxMC->GetX();
2589  Double_t mcProdVtxY = primVtxMC->GetY();
2590  Double_t mcProdVtxZ = primVtxMC->GetZ();
2591 
2592  // Process True Mesons
2593  if(TrueGammaCandidate0->GetV0Index()<fInputEvent->GetNumberOfV0s()){
2594  Bool_t isTrueMeson = kFALSE;
2595  Int_t gamma0MCLabel = -1;
2596  Int_t gamma0MotherLabel = -1;
2597  if (TrueGammaCandidate0->IsTrueConvertedPhoton()){
2598  gamma0MCLabel = TrueGammaCandidate0->GetMCParticleLabel(fMCEvent);
2599  if(gamma0MCLabel>-1){
2600  TParticle * gammaMC0 = (TParticle*)fMCEvent->Particle(gamma0MCLabel);
2601  gamma0MotherLabel=gammaMC0->GetFirstMother();
2602  }
2603  }
2604  if (!TrueGammaCandidate1->GetIsCaloPhoton()) AliFatal("CaloPhotonFlag has not been set. Aborting");
2605 
2606  Int_t gamma1MCLabel = TrueGammaCandidate1->GetCaloPhotonMCLabel(0); // get most probable MC label
2607  Int_t gamma1MotherLabel = -1;
2608  // check if
2609 
2610  TParticle * gammaMC1 = 0x0;
2611  if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
2612  // Daughters Gamma 1
2613  gammaMC1 = (TParticle*)fMCEvent->Particle(gamma1MCLabel);
2614  if (TrueGammaCandidate1->IsLargestComponentPhoton() || TrueGammaCandidate1->IsLargestComponentElectron()){ // largest component is electro magnetic
2615  // get mother of interest (pi0 or eta)
2616  if (TrueGammaCandidate1->IsLargestComponentPhoton()){ // for photons its the direct mother
2617  gamma1MotherLabel=gammaMC1->GetMother(0);
2618  }else if (TrueGammaCandidate1->IsLargestComponentElectron()){ // for electrons its either the direct mother or for conversions the grandmother
2619  if (TrueGammaCandidate1->IsConversion() && gammaMC1->GetMother(0)>-1) gamma1MotherLabel=fMCEvent->Particle(gammaMC1->GetMother(0))->GetMother(0);
2620  else gamma1MotherLabel=gammaMC1->GetMother(0);
2621  }
2622  }
2623  }
2624 
2625  if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
2626  if(((TParticle*)fMCEvent->Particle(gamma1MotherLabel))->GetPdgCode() == fMesonPDG){
2627  isTrueMeson=kTRUE;
2628  }
2629  }
2630 
2631  if(isTrueMeson ){// True Pion or Eta
2632  if (!matched){
2633  if (fHistoTrueMesonInvMassPt[fiCut]) fHistoTrueMesonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),fWeightJetJetMC);
2634  }else{
2635  if (fHistoTrueMesonMatchedInvMassPt[fiCut]) fHistoTrueMesonMatchedInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),fWeightJetJetMC);
2636  }
2637  if (fDoMesonQA > 0 && fIsMC < 2){
2638  if (TrueGammaCandidate1->IsLargestComponentPhoton() && !matched){
2639  if (fHistoTrueMesonCaloPhotonInvMassPt[fiCut]) fHistoTrueMesonCaloPhotonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2640  }
2641  if (TrueGammaCandidate1->IsLargestComponentElectron() && !matched){
2642  if (fHistoTrueMesonCaloElectronInvMassPt[fiCut]) fHistoTrueMesonCaloElectronInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2643  }
2644  if (TrueGammaCandidate1->IsLargestComponentElectron() && TrueGammaCandidate1->IsConversion() ){
2645  if ( !matched){
2646  fHistoTrueMesonCaloConvertedPhotonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2647  }
2648 
2649  if ((TrueGammaCandidate0->GetMCLabelPositive() == gamma1MCLabel || TrueGammaCandidate0->GetMCLabelNegative() == gamma1MCLabel) && isTrueMeson){
2650  if (fHistoTrueMesonCaloConvertedPhotonMatchedInvMassPt[fiCut])fHistoTrueMesonCaloConvertedPhotonMatchedInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2651  }
2652  }
2653  if ((TrueGammaCandidate1->IsMerged() || TrueGammaCandidate1->IsMergedPartConv() || TrueGammaCandidate1->IsDalitzMerged()) && !matched ){
2654  if (fHistoTrueMesonCaloMergedClusterInvMassPt[fiCut]) fHistoTrueMesonCaloMergedClusterInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2655  }
2656  if (TrueGammaCandidate1->IsMergedPartConv() && !matched){
2657  if (fHistoTrueMesonCaloMergedClusterPartConvInvMassPt[fiCut]) fHistoTrueMesonCaloMergedClusterPartConvInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2658  }
2659  }
2660  if (!matched){
2661  if (fDoMesonQA > 0){
2662  if (isTrueMeson){
2663  if (mesonCand->M() > fMesonInvMassWindow[0] && mesonCand->M() < fMesonInvMassWindow[1] ){
2664  if (fHistoTrueMesonPtY[fiCut]) fHistoTrueMesonPtY[fiCut]->Fill(mesonCand->Pt(),mesonCand->Rapidity()-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift(),fWeightJetJetMC);
2665  if (fHistoTrueMesonPtAlpha[fiCut]) fHistoTrueMesonPtAlpha[fiCut]->Fill(mesonCand->Pt(),mesonCand->GetAlpha(),fWeightJetJetMC);
2666  if (fHistoTrueMesonPtOpenAngle[fiCut]) fHistoTrueMesonPtOpenAngle[fiCut]->Fill(mesonCand->Pt(),mesonCand->GetOpeningAngle(),fWeightJetJetMC);
2667  if (fHistoTrueMotherMesonConvPhotonEtaPhi[fiCut])fHistoTrueMotherMesonConvPhotonEtaPhi[fiCut]->Fill(TrueGammaCandidate0->GetPhotonPhi(), TrueGammaCandidate0->GetPhotonEta(),fWeightJetJetMC);
2668  }
2669  }
2670  }
2671  Bool_t isPrimary = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsConversionPrimaryESD( fMCEvent, gamma0MotherLabel, mcProdVtxX, mcProdVtxY, mcProdVtxZ);
2672  if (isPrimary) {
2673  Float_t weighted= 1;
2674  if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(gamma1MotherLabel, fMCEvent, fInputEvent)){
2675  if (((TParticle*)fMCEvent->Particle(gamma1MotherLabel))->Pt()>0.005){
2676  weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(gamma1MotherLabel, fMCEvent, fInputEvent);
2677  // cout << "rec \t " <<gamma1MotherLabel << "\t" << weighted << endl;
2678  }
2679  }
2680  if (isTrueMeson){
2681  if (fHistoTruePrimaryMesonInvMassPt[fiCut]) fHistoTruePrimaryMesonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),weighted*fWeightJetJetMC);
2682  if (fHistoTruePrimaryMesonW0WeightingInvMassPt[fiCut]) fHistoTruePrimaryMesonW0WeightingInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),fWeightJetJetMC);
2683  if (fProfileTruePrimaryMesonWeightsInvMassPt[fiCut]) fProfileTruePrimaryMesonWeightsInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),weighted*fWeightJetJetMC);
2684  if (CheckVectorForDoubleCount(fVectorDoubleCountTrueMesons,gamma0MotherLabel)){
2685  if (fHistoDoubleCountTrueMesonInvMassPt[fiCut]) fHistoDoubleCountTrueMesonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),weighted*fWeightJetJetMC);
2686  FillMultipleCountMap(fMapMultipleCountTrueMesons,gamma0MotherLabel);
2687  }
2688  }
2689 
2690  if (fDoMesonQA > 0 && fIsMC < 2){
2691  if(isTrueMeson){ // Only primary pi0 for resolution
2692  if (fHistoTruePrimaryMesonMCPtResolPt[fiCut]) fHistoTruePrimaryMesonMCPtResolPt[fiCut]->Fill(((TParticle*)fMCEvent->Particle(gamma1MotherLabel))->Pt(),(mesonCand->Pt()-((TParticle*)fMCEvent->Particle(gamma1MotherLabel))->Pt())/((TParticle*)fMCEvent->Particle(gamma1MotherLabel))->Pt(),weighted*fWeightJetJetMC);
2693  }
2694  }
2695  }
2696  }
2697  }else if(!isTrueMeson ){ // Background
2698  if (fDoMesonQA > 1){
2699  if(gamma0MotherLabel>-1 && gamma1MotherLabel>-1){ // Both Tracks are Photons and have a mother but not Meson or Eta
2700  if (fHistoTrueBckGGInvMassPt[fiCut])fHistoTrueBckGGInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
2701  if(
2702  ( ((TParticle*)fMCEvent->Particle(gamma1MotherLabel))->GetPdgCode() == fMesonPDG
2703  && (TrueGammaCandidate1->IsMerged() || TrueGammaCandidate1->IsMergedPartConv()))
2704  ){
2705  if (fHistoTrueBckFullMesonContainedInOneClusterInvMassPt[fiCut]) fHistoTrueBckFullMesonContainedInOneClusterInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
2706  }else if( TrueGammaCandidate1->E()/mesonCand->E() > 0.7 ){
2707  if (fHistoTrueBckAsymEClustersInvMassPt[fiCut]) fHistoTrueBckAsymEClustersInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
2708  }
2709  }else { // No photon or without mother
2710  if (fHistoTrueBckContInvMassPt[fiCut]) fHistoTrueBckContInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
2711  }
2712  }
2713  }
2714  if (isTrueMeson && !matched){
2715  fVectorRecTrueMesons.push_back(gamma0MotherLabel);
2716  }
2717  }
2718 }
2719 
2720 //______________________________________________________________________
2721 void AliAnalysisTaskHeavyNeutralMesonToGG::ProcessTrueMesonCandidatesConvCaloAOD( AliAODConversionMother *mesonCand,
2722  AliAODConversionPhoton *TrueGammaCandidate0,
2723  AliAODConversionPhoton *TrueGammaCandidate1,
2724  Bool_t matched )
2725 {
2726  const AliVVertex* primVtxMC = fMCEvent->GetPrimaryVertex();
2727  Double_t mcProdVtxX = primVtxMC->GetX();
2728  Double_t mcProdVtxY = primVtxMC->GetY();
2729  Double_t mcProdVtxZ = primVtxMC->GetZ();
2730 
2731  // Process True Mesons
2732  TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
2733  if (AODMCTrackArray == NULL) return;
2734  Bool_t isTrueMeson = kFALSE;
2735 
2736  AliAODMCParticle *positiveMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate0->GetMCLabelPositive()));
2737  AliAODMCParticle *negativeMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate0->GetMCLabelNegative()));
2738 
2739  Int_t gamma0MCLabel = -1;
2740  Int_t gamma0MotherLabel = -1;
2741  if(!positiveMC||!negativeMC)
2742  return;
2743 
2744  if (TrueGammaCandidate0->IsTrueConvertedPhoton()){
2745  gamma0MCLabel = positiveMC->GetMother();
2746  AliAODMCParticle * gammaMC0 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma0MCLabel));
2747  gamma0MotherLabel=gammaMC0->GetMother();
2748  }
2749 
2750  if (!TrueGammaCandidate1->GetIsCaloPhoton()) AliFatal("CaloPhotonFlag has not been set. Aborting");
2751  Int_t gamma1MCLabel = TrueGammaCandidate1->GetCaloPhotonMCLabel(0); // get most probable MC label
2752  Int_t gamma1MotherLabel = -1;
2753  // check if
2754 
2755  AliAODMCParticle * gammaMC1 = 0x0;
2756  if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
2757  // Daughters Gamma 1
2758  gammaMC1 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MCLabel));
2759  if (TrueGammaCandidate1->IsLargestComponentPhoton() || TrueGammaCandidate1->IsLargestComponentElectron()){ // largest component is electro magnetic
2760  // get mother of interest (pi0 or eta)
2761  if (TrueGammaCandidate1->IsLargestComponentPhoton()){ // for photons its the direct mother
2762  gamma1MotherLabel=gammaMC1->GetMother();
2763  }else if (TrueGammaCandidate1->IsLargestComponentElectron()){ // for electrons its either the direct mother or for conversions the grandmother
2764  if (TrueGammaCandidate1->IsConversion()){
2765  AliAODMCParticle * gammaGrandMotherMC1 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gammaMC1->GetMother()));
2766  gamma1MotherLabel=gammaGrandMotherMC1->GetMother();
2767  }else gamma1MotherLabel=gammaMC1->GetMother();
2768  }
2769  }
2770  }
2771 
2772  if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
2773  if(((AliAODMCParticle*)AODMCTrackArray->At(gamma1MotherLabel))->GetPdgCode() == fMesonPDG){
2774  isTrueMeson=kTRUE;
2775  }
2776  }
2777 
2778  if(isTrueMeson ){// True Pion or Eta
2779  if (!matched){
2780  if (fHistoTrueMesonInvMassPt[fiCut]) fHistoTrueMesonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),fWeightJetJetMC);
2781 
2782  }else{
2783  if (fHistoTrueMesonMatchedInvMassPt[fiCut]) fHistoTrueMesonMatchedInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),fWeightJetJetMC);
2784  }
2785  if (fDoMesonQA > 0 && fIsMC < 2){
2786  if (TrueGammaCandidate1->IsLargestComponentPhoton() && !matched){
2787  if (fHistoTrueMesonCaloPhotonInvMassPt[fiCut]) fHistoTrueMesonCaloPhotonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2788  }
2789  if (TrueGammaCandidate1->IsLargestComponentElectron() && !matched) {
2790  if (fHistoTrueMesonCaloElectronInvMassPt[fiCut]) fHistoTrueMesonCaloElectronInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2791  }
2792  if (TrueGammaCandidate1->IsLargestComponentElectron() && TrueGammaCandidate1->IsConversion()){
2793  if (!matched)fHistoTrueMesonCaloConvertedPhotonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2794  if ((TrueGammaCandidate0->GetMCLabelPositive() == gamma1MCLabel || TrueGammaCandidate0->GetMCLabelNegative() == gamma1MCLabel) && isTrueMeson)
2795  if (fHistoTrueMesonCaloConvertedPhotonMatchedInvMassPt[fiCut])fHistoTrueMesonCaloConvertedPhotonMatchedInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2796  }
2797  if ((TrueGammaCandidate1->IsMerged() || TrueGammaCandidate1->IsMergedPartConv() || TrueGammaCandidate1->IsDalitzMerged()) && !matched ){
2798  if (fHistoTrueMesonCaloMergedClusterInvMassPt[fiCut]) fHistoTrueMesonCaloMergedClusterInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2799  }
2800  if (TrueGammaCandidate1->IsMergedPartConv() && !matched) {
2801  if (fHistoTrueMesonCaloMergedClusterPartConvInvMassPt[fiCut]) fHistoTrueMesonCaloMergedClusterPartConvInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2802  }
2803  }
2804 
2805  if ( !matched){
2806  if (fDoMesonQA > 0){
2807  if (mesonCand->M() > fMesonInvMassWindow[0] && mesonCand->M() < fMesonInvMassWindow[1] ){
2808  if (fHistoTrueMesonPtY[fiCut]) fHistoTrueMesonPtY[fiCut]->Fill(mesonCand->Pt(),mesonCand->Rapidity()-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift(),fWeightJetJetMC);
2809  if (fHistoTrueMesonPtAlpha[fiCut]) fHistoTrueMesonPtAlpha[fiCut]->Fill(mesonCand->Pt(),mesonCand->GetAlpha(),fWeightJetJetMC);
2810  if (fHistoTrueMesonPtOpenAngle[fiCut]) fHistoTrueMesonPtOpenAngle[fiCut]->Fill(mesonCand->Pt(),mesonCand->GetOpeningAngle(),fWeightJetJetMC);
2811  if (fHistoTrueMotherMesonConvPhotonEtaPhi[fiCut]) fHistoTrueMotherMesonConvPhotonEtaPhi[fiCut]->Fill(TrueGammaCandidate0->GetPhotonPhi(), TrueGammaCandidate0->GetPhotonEta(),fWeightJetJetMC);
2812  }
2813  }
2814 
2815  Bool_t isPrimary = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsConversionPrimaryAOD(fInputEvent, static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma0MotherLabel)), mcProdVtxX, mcProdVtxY, mcProdVtxZ);
2816  if(isPrimary){
2817  // Only primary pi0 for efficiency calculation
2818  Float_t weighted= 1;
2819  if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(gamma1MotherLabel, 0x0, fInputEvent)){
2820  if (static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt()>0.005){
2821  weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(gamma1MotherLabel, 0x0, fInputEvent);
2822  // cout << "rec \t " <<gamma1MotherLabel << "\t" << weighted << endl;
2823  }
2824  }
2825  if (isTrueMeson){
2826  if (fHistoTruePrimaryMesonInvMassPt[fiCut]) fHistoTruePrimaryMesonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),weighted*fWeightJetJetMC);
2827  if (fHistoTruePrimaryMesonW0WeightingInvMassPt[fiCut]) fHistoTruePrimaryMesonW0WeightingInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),fWeightJetJetMC);
2828  if (fProfileTruePrimaryMesonWeightsInvMassPt[fiCut]) fProfileTruePrimaryMesonWeightsInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),weighted*fWeightJetJetMC);
2829  if (CheckVectorForDoubleCount(fVectorDoubleCountTrueMesons,gamma0MotherLabel)){
2830  if (fHistoDoubleCountTrueMesonInvMassPt[fiCut]) fHistoDoubleCountTrueMesonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),weighted*fWeightJetJetMC);
2831  FillMultipleCountMap(fMapMultipleCountTrueMesons,gamma0MotherLabel);
2832  }
2833  }
2834  if (fDoMesonQA > 0 && fIsMC < 2){
2835  if(isTrueMeson){ // Only primary pi0 for resolution
2836  if (fHistoTruePrimaryMesonMCPtResolPt[fiCut]) fHistoTruePrimaryMesonMCPtResolPt[fiCut]->Fill(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),
2837  (mesonCand->Pt()-static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt())/static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),weighted*fWeightJetJetMC);
2838 
2839  }
2840  }
2841  }
2842  }
2843  }else if(!isTrueMeson ) { // Background
2844  if (fDoMesonQA > 1){
2845  if(gamma0MotherLabel>-1 && gamma1MotherLabel>-1){ // Both Tracks are Photons and have a mother but not Meson or Eta
2846  if (fHistoTrueBckGGInvMassPt[fiCut]) fHistoTrueBckGGInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2847 
2848  if( (((AliAODMCParticle*)AODMCTrackArray->At(gamma1MotherLabel))->GetPdgCode() == fMesonPDG
2849  && ((TrueGammaCandidate1->IsMerged() || TrueGammaCandidate1->IsMergedPartConv())))
2850  ){
2851  if (fHistoTrueBckFullMesonContainedInOneClusterInvMassPt[fiCut]) fHistoTrueBckFullMesonContainedInOneClusterInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
2852  }else if( TrueGammaCandidate1->E()/mesonCand->E() > 0.7 ){
2853  if (fHistoTrueBckAsymEClustersInvMassPt[fiCut]) fHistoTrueBckAsymEClustersInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
2854  }
2855  }else { // No photon or without mother
2856  if (fHistoTrueBckContInvMassPt[fiCut]) fHistoTrueBckContInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2857  }
2858  }
2859  }
2860  if (isTrueMeson && !matched){
2861  fVectorRecTrueMesons.push_back(gamma0MotherLabel);
2862  }
2863 }
2864 
2865 
2866 //______________________________________________________________________
2867 void AliAnalysisTaskHeavyNeutralMesonToGG::ProcessTrueMesonCandidatesCalo( AliAODConversionMother *mesonCand,
2868  AliAODConversionPhoton *TrueGammaCandidate0,
2869  AliAODConversionPhoton *TrueGammaCandidate1)
2870 {
2871  // Process True Mesons
2872  const AliVVertex* primVtxMC = fMCEvent->GetPrimaryVertex();
2873  Double_t mcProdVtxX = primVtxMC->GetX();
2874  Double_t mcProdVtxY = primVtxMC->GetY();
2875  Double_t mcProdVtxZ = primVtxMC->GetZ();
2876 
2877  Bool_t isTrueMeson = kFALSE;
2878  Int_t gamma0MCLabel = TrueGammaCandidate0->GetCaloPhotonMCLabel(0); // get most probable MC label
2879  Int_t gamma0MotherLabel = -1;
2880 
2881  TParticle * gammaMC0 = 0x0;
2882  if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
2883  gammaMC0 = (TParticle*)fMCEvent->Particle(gamma0MCLabel);
2884  if (TrueGammaCandidate0->IsLargestComponentPhoton() || TrueGammaCandidate0->IsLargestComponentElectron()){ // largest component is electro magnetic
2885  // get mother of interest (pi0 or eta)
2886  if (TrueGammaCandidate0->IsLargestComponentPhoton()){ // for photons its the direct mother
2887  gamma0MotherLabel=gammaMC0->GetMother(0);
2888  } else if (TrueGammaCandidate0->IsLargestComponentElectron()){ // for electrons its either the direct mother or for conversions the grandmother
2889  if (TrueGammaCandidate0->IsConversion() && (gammaMC0->GetMother(0) > -1)){
2890  gamma0MotherLabel=fMCEvent->Particle(gammaMC0->GetMother(0))->GetMother(0);
2891  } else {
2892  gamma0MotherLabel=gammaMC0->GetMother(0);
2893  }
2894  }
2895  }
2896  }
2897  if (!TrueGammaCandidate1->GetIsCaloPhoton()) AliFatal("CaloPhotonFlag has not been set. Aborting");
2898 
2899  Int_t gamma1MCLabel = TrueGammaCandidate1->GetCaloPhotonMCLabel(0); // get most probable MC label
2900  Int_t gamma1MotherLabel = -1;
2901  // check if
2902 
2903  TParticle * gammaMC1 = 0x0;
2904  if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
2905  // Daughters Gamma 1
2906  gammaMC1 = (TParticle*)fMCEvent->Particle(gamma1MCLabel);
2907  if (TrueGammaCandidate1->IsLargestComponentPhoton() || TrueGammaCandidate1->IsLargestComponentElectron()){ // largest component is electro magnetic
2908  // get mother of interest (pi0 or eta)
2909  if (TrueGammaCandidate1->IsLargestComponentPhoton()){ // for photons its the direct mother
2910  gamma1MotherLabel = gammaMC1->GetMother(0);
2911  } else if (TrueGammaCandidate1->IsLargestComponentElectron()){ // for electrons its either the direct mother or for conversions the grandmother
2912  if (TrueGammaCandidate1->IsConversion()){
2913  if(gammaMC1->GetMother(0) > -1) gamma1MotherLabel = fMCEvent->Particle(gammaMC1->GetMother(0))->GetMother(0);
2914  } else {
2915  gamma1MotherLabel=gammaMC1->GetMother(0);
2916  }
2917  }
2918  }
2919  }
2920 
2921  if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
2922  if(((TParticle*)fMCEvent->Particle(gamma1MotherLabel))->GetPdgCode() == fMesonPDG){
2923  isTrueMeson=kTRUE;
2924  }
2925  }
2926 
2927  if(isTrueMeson ){// True Meson
2928  if (fHistoTrueMesonInvMassPt[fiCut]) fHistoTrueMesonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
2929  if (fDoMesonQA > 0 && fDoMesonQA < 3 && fIsMC < 2){
2930  // both gammas are real gammas
2931  if (TrueGammaCandidate0->IsLargestComponentPhoton() && TrueGammaCandidate1->IsLargestComponentPhoton()) {
2932  if (fHistoTrueMesonCaloPhotonInvMassPt[fiCut]) fHistoTrueMesonCaloPhotonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2933  }
2934  // both particles are electrons
2935  if (TrueGammaCandidate0->IsLargestComponentElectron() && TrueGammaCandidate1->IsLargestComponentElectron() ) {
2936  if (fHistoTrueMesonCaloElectronInvMassPt[fiCut]) fHistoTrueMesonCaloElectronInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2937  }
2938  // both particles are converted electrons
2939  if ((TrueGammaCandidate0->IsLargestComponentElectron() && TrueGammaCandidate0->IsConversion()) && (TrueGammaCandidate1->IsLargestComponentElectron() && TrueGammaCandidate1->IsConversion()) ){
2940  fHistoTrueMesonCaloConvertedPhotonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2941  }
2942  // 1 gamma is converted the other one is normals
2943  if ( (TrueGammaCandidate0->IsLargestComponentPhoton() && TrueGammaCandidate1->IsLargestComponentElectron() && TrueGammaCandidate1->IsConversion()) ||
2944  (TrueGammaCandidate1->IsLargestComponentPhoton() && TrueGammaCandidate0->IsLargestComponentElectron() && TrueGammaCandidate0->IsConversion())
2945  ) {
2946  fHistoTrueMesonCaloMixedPhotonConvPhotonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2947  }
2948 
2949  // at least one of the photon is merged
2950  if (TrueGammaCandidate0->IsMerged() || TrueGammaCandidate0->IsMergedPartConv() || TrueGammaCandidate0->IsDalitzMerged() || TrueGammaCandidate1->IsMerged() || TrueGammaCandidate1->IsMergedPartConv() || TrueGammaCandidate1->IsDalitzMerged() ){
2951  if (fHistoTrueMesonCaloMergedClusterInvMassPt[fiCut]) fHistoTrueMesonCaloMergedClusterInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2952  }
2953  // at least one of the photon is merged and part conv
2954  if (TrueGammaCandidate1->IsMergedPartConv() || TrueGammaCandidate0->IsMergedPartConv()) {
2955  if (fHistoTrueMesonCaloMergedClusterPartConvInvMassPt[fiCut]) fHistoTrueMesonCaloMergedClusterPartConvInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
2956  }
2957  }
2958 
2959  if (fDoMesonQA > 0 && fDoMesonQA < 3){
2960  if ( mesonCand->M() > fMesonInvMassWindow[0] && mesonCand->M() < fMesonInvMassWindow[1] ){
2961  if (fHistoTrueMesonPtY[fiCut]) fHistoTrueMesonPtY[fiCut]->Fill(mesonCand->Pt(),mesonCand->Rapidity()-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift(), fWeightJetJetMC);
2962  if (fHistoTrueMesonPtAlpha[fiCut]) fHistoTrueMesonPtAlpha[fiCut]->Fill(mesonCand->Pt(),TMath::Abs(mesonCand->GetAlpha()), fWeightJetJetMC);
2963  if (fHistoTrueMesonPtOpenAngle[fiCut]) fHistoTrueMesonPtOpenAngle[fiCut]->Fill(mesonCand->Pt(),mesonCand->GetOpeningAngle(), fWeightJetJetMC);
2964  }
2965 
2966  }
2967  Bool_t isPrimary = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsConversionPrimaryESD( fMCEvent, gamma0MotherLabel, mcProdVtxX, mcProdVtxY, mcProdVtxZ);
2968  if(isPrimary){ // Only primary pi0 for efficiency calculation
2969  // filling primary histograms
2970  Float_t weighted= 1;
2971  if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(gamma1MotherLabel, fMCEvent, fInputEvent)){
2972  if (((TParticle*)fMCEvent->Particle(gamma1MotherLabel))->Pt()>0.005){
2973  weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(gamma1MotherLabel, fMCEvent, fInputEvent);
2974  // cout << "rec \t " <<gamma1MotherLabel << "\t" << weighted << endl;
2975  }
2976  }
2977  if (isTrueMeson){
2978  if (fHistoTruePrimaryMesonInvMassPt[fiCut]) fHistoTruePrimaryMesonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),weighted* fWeightJetJetMC);
2979  if (fHistoTruePrimaryMesonW0WeightingInvMassPt[fiCut]) fHistoTruePrimaryMesonW0WeightingInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
2980  if (fProfileTruePrimaryMesonWeightsInvMassPt[fiCut]) fProfileTruePrimaryMesonWeightsInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),weighted* fWeightJetJetMC);
2981  if (CheckVectorForDoubleCount(fVectorDoubleCountTrueMesons,gamma0MotherLabel) && fHistoDoubleCountTrueMesonInvMassPt[fiCut]) fHistoDoubleCountTrueMesonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), weighted*fWeightJetJetMC);
2982  }
2983  if (fDoMesonQA > 0 && fDoMesonQA < 3 && fIsMC<2){
2984  if(isTrueMeson){ // Only primary pi0 for resolution
2985  if (fHistoTruePrimaryMesonMCPtResolPt[fiCut]) fHistoTruePrimaryMesonMCPtResolPt[fiCut]->Fill(((TParticle*)fMCEvent->Particle(gamma1MotherLabel))->Pt(),(mesonCand->Pt()-((TParticle*)fMCEvent->Particle(gamma1MotherLabel))->Pt())/((TParticle*)fMCEvent->Particle(gamma1MotherLabel))->Pt(),weighted* fWeightJetJetMC);
2986  }
2987  }
2988  }
2989  } else if(!isTrueMeson ){ // Background
2990  if (fDoMesonQA > 1 && fDoMesonQA < 3){
2991  if(gamma0MotherLabel>-1 && gamma1MotherLabel>-1){ // Both Tracks are Photons and have a mother but not Pi0 or Eta
2992  if (fHistoTrueBckGGInvMassPt[fiCut]) fHistoTrueBckGGInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
2993 
2994  if( (((TParticle*)fMCEvent->Particle(gamma0MotherLabel))->GetPdgCode() == fMesonPDG
2995  && (TrueGammaCandidate0->IsMerged() || TrueGammaCandidate0->IsMergedPartConv()))
2996  ||
2997  (((TParticle*)fMCEvent->Particle(gamma1MotherLabel))->GetPdgCode() == fMesonPDG
2998  && (TrueGammaCandidate1->IsMerged() || TrueGammaCandidate1->IsMergedPartConv()))
2999  ){
3000  if (fHistoTrueBckFullMesonContainedInOneClusterInvMassPt[fiCut]) fHistoTrueBckFullMesonContainedInOneClusterInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
3001  }else if( (TrueGammaCandidate0->E()/mesonCand->E() > 0.7) || (TrueGammaCandidate1->E()/mesonCand->E() > 0.7) ){
3002  if (fHistoTrueBckAsymEClustersInvMassPt[fiCut]) fHistoTrueBckAsymEClustersInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
3003  }
3004  } else { // No photon or without mother
3005  if (fHistoTrueBckContInvMassPt[fiCut]) fHistoTrueBckContInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
3006  }
3007  }
3008  }
3009 }
3010 
3011 //______________________________________________________________________
3012 void AliAnalysisTaskHeavyNeutralMesonToGG::ProcessTrueMesonCandidatesCaloAOD( AliAODConversionMother *mesonCand,
3013  AliAODConversionPhoton *TrueGammaCandidate0,
3014  AliAODConversionPhoton *TrueGammaCandidate1)
3015 {
3016  const AliVVertex* primVtxMC = fMCEvent->GetPrimaryVertex();
3017  Double_t mcProdVtxX = primVtxMC->GetX();
3018  Double_t mcProdVtxY = primVtxMC->GetY();
3019  Double_t mcProdVtxZ = primVtxMC->GetZ();
3020 
3021  // Process True Mesons
3022  TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
3023  if (AODMCTrackArray == NULL) return;
3024 
3025  Bool_t isTrueMeson = kFALSE;
3026  Int_t gamma0MCLabel = TrueGammaCandidate0->GetCaloPhotonMCLabel(0); // get most probable MC label
3027  Int_t gamma0MotherLabel = -1;
3028 
3029  // check if
3030  AliAODMCParticle * gammaMC0 = 0x0;
3031  if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
3032  // Daughters Gamma 0
3033  gammaMC0 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma0MCLabel));
3034  if (TrueGammaCandidate0->IsLargestComponentPhoton() || TrueGammaCandidate0->IsLargestComponentElectron()){ // largest component is electro magnetic
3035  // get mother of interest (pi0 or eta)
3036  if (TrueGammaCandidate0->IsLargestComponentPhoton()){ // for photons its the direct mother
3037  gamma0MotherLabel=gammaMC0->GetMother();
3038  } else if (TrueGammaCandidate0->IsLargestComponentElectron()){ // for electrons its either the direct mother or for conversions the grandmother
3039  if (TrueGammaCandidate0->IsConversion()){
3040  AliAODMCParticle * gammaGrandMotherMC0 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gammaMC0->GetMother()));
3041  gamma0MotherLabel=gammaGrandMotherMC0->GetMother();
3042  } else gamma0MotherLabel=gammaMC0->GetMother();
3043  }
3044  }
3045  }
3046 
3047  Int_t gamma1MCLabel = TrueGammaCandidate1->GetCaloPhotonMCLabel(0); // get most probable MC label
3048  Int_t gamma1MotherLabel = -1;
3049 
3050  // check if
3051  AliAODMCParticle *gammaMC1 = 0x0;
3052  if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
3053  // Daughters Gamma 1
3054  gammaMC1 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MCLabel));
3055  if (TrueGammaCandidate1->IsLargestComponentPhoton() || TrueGammaCandidate1->IsLargestComponentElectron()){ // largest component is electro magnetic
3056  // get mother of interest (pi0 or eta)
3057  if (TrueGammaCandidate1->IsLargestComponentPhoton()){ // for photons its the direct mother
3058  gamma1MotherLabel=gammaMC1->GetMother();
3059  } else if (TrueGammaCandidate1->IsLargestComponentElectron()){ // for electrons its either the direct mother or for conversions the grandmother
3060  if (TrueGammaCandidate1->IsConversion()){
3061  AliAODMCParticle * gammaGrandMotherMC1 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gammaMC1->GetMother()));
3062  gamma1MotherLabel=gammaGrandMotherMC1->GetMother();
3063  } else gamma1MotherLabel=gammaMC1->GetMother();
3064  }
3065  }
3066  }
3067 
3068  if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
3069  if(((AliAODMCParticle*)AODMCTrackArray->At(gamma1MotherLabel))->GetPdgCode() == fMesonPDG){
3070  isTrueMeson=kTRUE;
3071  }
3072  }
3073 
3074  if(isTrueMeson ){// True Meson
3075  if (fHistoTrueMesonInvMassPt[fiCut]) fHistoTrueMesonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
3076  if (fDoMesonQA > 0 && fDoMesonQA < 3 && fIsMC < 2){
3077  // both gammas are real gammas
3078  if (TrueGammaCandidate0->IsLargestComponentPhoton() && TrueGammaCandidate1->IsLargestComponentPhoton()) {
3079  if (fHistoTrueMesonCaloPhotonInvMassPt[fiCut]) fHistoTrueMesonCaloPhotonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
3080  }
3081  // both particles are electrons
3082  if (TrueGammaCandidate0->IsLargestComponentElectron() && TrueGammaCandidate1->IsLargestComponentElectron() ) {
3083  if (fHistoTrueMesonCaloElectronInvMassPt[fiCut]) fHistoTrueMesonCaloElectronInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
3084  }
3085  // both particles are converted electrons
3086  if ((TrueGammaCandidate0->IsLargestComponentElectron() && TrueGammaCandidate0->IsConversion()) && (TrueGammaCandidate1->IsLargestComponentElectron() && TrueGammaCandidate1->IsConversion()) ){
3087  fHistoTrueMesonCaloConvertedPhotonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
3088  }
3089  // 1 gamma is converted the other one is normals
3090  if ( (TrueGammaCandidate0->IsLargestComponentPhoton() && TrueGammaCandidate1->IsLargestComponentElectron() && TrueGammaCandidate1->IsConversion()) ||
3091  (TrueGammaCandidate1->IsLargestComponentPhoton() && TrueGammaCandidate0->IsLargestComponentElectron() && TrueGammaCandidate0->IsConversion())
3092  ) {
3093  fHistoTrueMesonCaloMixedPhotonConvPhotonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
3094  }
3095 
3096  // at least one of the photon is merged
3097  if (TrueGammaCandidate0->IsMerged() || TrueGammaCandidate0->IsMergedPartConv() || TrueGammaCandidate0->IsDalitzMerged() || TrueGammaCandidate1->IsMerged() || TrueGammaCandidate1->IsMergedPartConv() || TrueGammaCandidate1->IsDalitzMerged() ){
3098  if (fHistoTrueMesonCaloMergedClusterInvMassPt[fiCut]) fHistoTrueMesonCaloMergedClusterInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
3099  }
3100  // at least one of the photon is merged and part conv
3101  if (TrueGammaCandidate1->IsMergedPartConv() || TrueGammaCandidate0->IsMergedPartConv()) {
3102  if (fHistoTrueMesonCaloMergedClusterPartConvInvMassPt[fiCut]) fHistoTrueMesonCaloMergedClusterPartConvInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
3103  }
3104  }
3105 
3106  if (fDoMesonQA > 0 && fDoMesonQA < 3){
3107  if ( mesonCand->M() > fMesonInvMassWindow[0] && mesonCand->M() < fMesonInvMassWindow[1] ){
3108  if (fHistoTrueMesonPtY[fiCut]) fHistoTrueMesonPtY[fiCut]->Fill(mesonCand->Pt(),mesonCand->Rapidity()-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift(), fWeightJetJetMC);
3109  if (fHistoTrueMesonPtAlpha[fiCut]) fHistoTrueMesonPtAlpha[fiCut]->Fill(mesonCand->Pt(),TMath::Abs(mesonCand->GetAlpha()), fWeightJetJetMC);
3110  if (fHistoTrueMesonPtOpenAngle[fiCut]) fHistoTrueMesonPtOpenAngle[fiCut]->Fill(mesonCand->Pt(),mesonCand->GetOpeningAngle(), fWeightJetJetMC);
3111  }
3112  }
3113 
3114  Bool_t isPrimary = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsConversionPrimaryAOD(fInputEvent, static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma0MotherLabel)), mcProdVtxX, mcProdVtxY, mcProdVtxZ);
3115  if(isPrimary){ // Only primary pi0 for efficiency calculation
3116  Float_t weighted= 1;
3117  if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(gamma1MotherLabel, 0x0, fInputEvent)){
3118  if (static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt()>0.005){
3119  weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(gamma1MotherLabel, 0x0, fInputEvent);
3120  }
3121  }
3122 
3123  if (fHistoTruePrimaryMesonInvMassPt[fiCut]) fHistoTruePrimaryMesonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),weighted* fWeightJetJetMC);
3124  if (fHistoTruePrimaryMesonW0WeightingInvMassPt[fiCut]) fHistoTruePrimaryMesonW0WeightingInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
3125  if (fProfileTruePrimaryMesonWeightsInvMassPt[fiCut]) fProfileTruePrimaryMesonWeightsInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),weighted* fWeightJetJetMC);
3126  if (CheckVectorForDoubleCount(fVectorDoubleCountTrueMesons,gamma0MotherLabel) && fHistoDoubleCountTrueMesonInvMassPt[fiCut]) fHistoDoubleCountTrueMesonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), weighted*fWeightJetJetMC);
3127  if (fDoMesonQA > 0 && fDoMesonQA < 3 && fIsMC<2){
3128  if (fHistoTruePrimaryMesonMCPtResolPt[fiCut]) fHistoTruePrimaryMesonMCPtResolPt[fiCut]->Fill(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),
3129  (mesonCand->Pt()-static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt())/static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),weighted* fWeightJetJetMC);
3130  }
3131  }
3132  } else if(!isTrueMeson ) { // Background
3133  if (fDoMesonQA > 1 && fDoMesonQA < 3){
3134  if(gamma0MotherLabel>-1 && gamma1MotherLabel>-1){ // Both Tracks are Photons and have a mother but not Pi0 or Eta
3135  if (fHistoTrueBckGGInvMassPt[fiCut]) fHistoTrueBckGGInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
3136 
3137  if( (((AliAODMCParticle*)AODMCTrackArray->At(gamma0MotherLabel))->GetPdgCode() == fMesonPDG
3138  && (TrueGammaCandidate0->IsMerged() || TrueGammaCandidate0->IsMergedPartConv()))
3139  ||
3140  (((AliAODMCParticle*)AODMCTrackArray->At(gamma1MotherLabel))->GetPdgCode() == fMesonPDG
3141  && (TrueGammaCandidate1->IsMerged() || TrueGammaCandidate1->IsMergedPartConv()))
3142  ){
3143  if (fHistoTrueBckFullMesonContainedInOneClusterInvMassPt[fiCut]) fHistoTrueBckFullMesonContainedInOneClusterInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
3144  }else if( (TrueGammaCandidate0->E()/mesonCand->E() > 0.7) || (TrueGammaCandidate1->E()/mesonCand->E() > 0.7) ){
3145  if (fHistoTrueBckAsymEClustersInvMassPt[fiCut]) fHistoTrueBckAsymEClustersInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
3146  }
3147  } else { // No photon or without mother
3148  if (fHistoTrueBckContInvMassPt[fiCut]) fHistoTrueBckContInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(), fWeightJetJetMC);
3149  }
3150  }
3151  }
3152 }
3153 
3154 
3155 //______________________________________________________________________
3156 void AliAnalysisTaskHeavyNeutralMesonToGG::ProcessTrueMesonCandidatesConv( AliAODConversionMother *mesonCand,
3157  AliAODConversionPhoton *TrueGammaCandidate0,
3158  AliAODConversionPhoton *TrueGammaCandidate1)
3159 {
3160  // Process True Mesons
3161  const AliVVertex* primVtxMC = fMCEvent->GetPrimaryVertex();
3162  Double_t mcProdVtxX = primVtxMC->GetX();
3163  Double_t mcProdVtxY = primVtxMC->GetY();
3164  Double_t mcProdVtxZ = primVtxMC->GetZ();
3165 
3166  if(TrueGammaCandidate0->GetV0Index()<fInputEvent->GetNumberOfV0s()){
3167  Bool_t isTrueMeson = kFALSE;
3168  Bool_t isTrueMesonDalitz = kFALSE;
3169  Bool_t gamma0DalitzCand = kFALSE;
3170  Bool_t gamma1DalitzCand = kFALSE;
3171  Int_t gamma0MCLabel = TrueGammaCandidate0->GetMCParticleLabel(fMCEvent);
3172  Int_t gamma0MotherLabel = -1;
3173  if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
3174  // Daughters Gamma 0
3175  TParticle * negativeMC = (TParticle*)TrueGammaCandidate0->GetNegativeMCDaughter(fMCEvent);
3176  TParticle * positiveMC = (TParticle*)TrueGammaCandidate0->GetPositiveMCDaughter(fMCEvent);
3177  TParticle * gammaMC0 = (TParticle*)fMCEvent->Particle(gamma0MCLabel);
3178  if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){ // Electrons ...
3179  if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
3180  if(gammaMC0->GetPdgCode() == 22){ // ... with Gamma Mother
3181  gamma0MotherLabel=gammaMC0->GetFirstMother();
3182  }
3183  }
3184  if(gammaMC0->GetPdgCode() ==fMesonPDG){ // Dalitz candidate
3185  gamma0DalitzCand = kTRUE;
3186  gamma0MotherLabel=-fMesonPDG;
3187  }
3188  }
3189  }
3190  if(TrueGammaCandidate1->GetV0Index()<fInputEvent->GetNumberOfV0s()){
3191  Int_t gamma1MCLabel = TrueGammaCandidate1->GetMCParticleLabel(fMCEvent);
3192  Int_t gamma1MotherLabel = -1;
3193  if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
3194  // Daughters Gamma 1
3195  TParticle * negativeMC = (TParticle*)TrueGammaCandidate1->GetNegativeMCDaughter(fMCEvent);
3196  TParticle * positiveMC = (TParticle*)TrueGammaCandidate1->GetPositiveMCDaughter(fMCEvent);
3197  TParticle * gammaMC1 = (TParticle*)fMCEvent->Particle(gamma1MCLabel);
3198  if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){ // Electrons ...
3199  if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
3200  if(gammaMC1->GetPdgCode() == 22){ // ... with Gamma Mother
3201  gamma1MotherLabel=gammaMC1->GetFirstMother();
3202  }
3203  }
3204  if(gammaMC1->GetPdgCode() ==fMesonPDG ){ // Dalitz candidate
3205  gamma1DalitzCand = kTRUE;
3206  gamma1MotherLabel=-fMesonPDG;
3207  }
3208  }
3209  }
3210  if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
3211  if(((TParticle*)fMCEvent->Particle(gamma1MotherLabel))->GetPdgCode() == fMesonPDG){
3212  isTrueMeson=kTRUE;
3213  if (CheckVectorForDoubleCount(fVectorDoubleCountTrueMesons,gamma0MotherLabel)){
3214  if (fHistoDoubleCountTrueMesonInvMassPt[fiCut]) fHistoDoubleCountTrueMesonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),fWeightJetJetMC);
3215  FillMultipleCountMap(fMapMultipleCountTrueMesons,gamma0MotherLabel);
3216  }
3217  }
3218  }
3219 
3220  //Identify Dalitz candidate
3221  if (gamma1DalitzCand || gamma0DalitzCand){
3222  if (gamma0DalitzCand && gamma0MCLabel >=0 && gamma0MCLabel==gamma1MotherLabel){
3223  if (gamma0MotherLabel == -fMesonPDG) isTrueMesonDalitz = kTRUE;
3224  }
3225  if (gamma1DalitzCand && gamma1MCLabel >=0 && gamma1MCLabel==gamma0MotherLabel){
3226  if (gamma1MotherLabel == -fMesonPDG) isTrueMesonDalitz = kTRUE;
3227  }
3228  }
3229 
3230 
3231  if(isTrueMeson ){// True Meosn
3232  if (fHistoTrueMesonInvMassPt[fiCut])fHistoTrueMesonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),fWeightJetJetMC);
3233  if (fDoMesonQA > 0){
3234  if ( mesonCand->M() > fMesonInvMassWindow[0] && mesonCand->M() < fMesonInvMassWindow[1]){
3235  if(fIsMC < 2){
3236  if (fHistoTrueMesonPtY[fiCut]) fHistoTrueMesonPtY[fiCut]->Fill(mesonCand->Pt(),mesonCand->Rapidity()-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift());
3237  if (fHistoTrueMesonPtOpenAngle[fiCut]) fHistoTrueMesonPtOpenAngle[fiCut]->Fill(mesonCand->Pt(),mesonCand->GetOpeningAngle());
3238  if (fHistoTrueMesonPtAlpha[fiCut]) fHistoTrueMesonPtAlpha[fiCut]->Fill(mesonCand->Pt(),TMath::Abs(mesonCand->GetAlpha()),fWeightJetJetMC);
3239  }
3240  }
3241  }
3242  Bool_t isPrimary = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsConversionPrimaryESD( fMCEvent, gamma0MotherLabel, mcProdVtxX, mcProdVtxY, mcProdVtxZ);
3243  if(isPrimary){ // Only primary pi0 for efficiency calculation
3244  Float_t weighted= 1;
3245  if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(gamma1MotherLabel, fMCEvent, fInputEvent)){
3246  if (((TParticle*)fMCEvent->Particle(gamma1MotherLabel))->Pt()>0.005){
3247  weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(gamma1MotherLabel, fMCEvent, fInputEvent);
3248  // cout << "rec \t " <<gamma1MotherLabel << "\t" << weighted << endl;
3249  }
3250  }
3251  if (fHistoTruePrimaryMesonInvMassPt[fiCut]) fHistoTruePrimaryMesonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),weighted*fWeightJetJetMC);
3252  if (fHistoTruePrimaryMesonW0WeightingInvMassPt[fiCut]) fHistoTruePrimaryMesonW0WeightingInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),fWeightJetJetMC);
3253  if (fProfileTruePrimaryMesonWeightsInvMassPt[fiCut]) fProfileTruePrimaryMesonWeightsInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),weighted*fWeightJetJetMC);
3254 
3255  if (fDoMesonQA > 0 && fIsMC < 2){
3256  if (fHistoTruePrimaryMesonMCPtResolPt[fiCut]) fHistoTruePrimaryMesonMCPtResolPt[fiCut]->Fill(((TParticle*)fMCEvent->Particle(gamma1MotherLabel))->Pt(),(mesonCand->Pt()-((TParticle*)fMCEvent->Particle(gamma1MotherLabel))->Pt())/((TParticle*)fMCEvent->Particle(gamma1MotherLabel))->Pt(),weighted);
3257  }
3258  }
3259  } else if(!isTrueMeson ){ // Background
3260  if (fDoMesonQA > 1 && fIsMC < 2){
3261  if(gamma0MotherLabel>-1 && gamma1MotherLabel>-1){ // Both Tracks are Photons and have a mother but not Pi0 or Eta
3262  if (fHistoTrueBckGGInvMassPt[fiCut])fHistoTrueBckGGInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
3263  } else { // No photon or without mother
3264  if (fHistoTrueBckContInvMassPt[fiCut]) fHistoTrueBckContInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
3265  }
3266  }
3267  if (!isTrueMesonDalitz && (gamma0DalitzCand || gamma1DalitzCand)){
3268  if (fDoMesonQA > 1 && fIsMC < 2 && fHistoTrueBckContInvMassPt[fiCut]) fHistoTrueBckContInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
3269  }
3270  }
3271  }
3272  }
3273 }
3274 
3275 //______________________________________________________________________
3276 void AliAnalysisTaskHeavyNeutralMesonToGG::ProcessTrueMesonCandidatesConvAOD( AliAODConversionMother *mesonCand,
3277  AliAODConversionPhoton *TrueGammaCandidate0,
3278  AliAODConversionPhoton *TrueGammaCandidate1)
3279 {
3280  const AliVVertex* primVtxMC = fMCEvent->GetPrimaryVertex();
3281  Double_t mcProdVtxX = primVtxMC->GetX();
3282  Double_t mcProdVtxY = primVtxMC->GetY();
3283  Double_t mcProdVtxZ = primVtxMC->GetZ();
3284 
3285  // Process True Mesons
3286  TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
3287  Bool_t isTrueMeson = kFALSE;
3288  Bool_t isTrueMesonDalitz = kFALSE;
3289  Bool_t gamma0DalitzCand = kFALSE;
3290  Bool_t gamma1DalitzCand = kFALSE;
3291 
3292  if (AODMCTrackArray!=NULL && TrueGammaCandidate0 != NULL){
3293  AliAODMCParticle *positiveMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate0->GetMCLabelPositive()));
3294  AliAODMCParticle *negativeMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate0->GetMCLabelNegative()));
3295 
3296  Int_t gamma0MCLabel = -1;
3297  Int_t gamma0MotherLabel = -1;
3298  if(!positiveMC||!negativeMC)
3299  return;
3300 
3301  if(positiveMC->GetMother()>-1&&(negativeMC->GetMother() == positiveMC->GetMother())){
3302  gamma0MCLabel = positiveMC->GetMother();
3303  }
3304 
3305  if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
3306  // Daughters Gamma 0
3307  AliAODMCParticle * gammaMC0 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma0MCLabel));
3308  if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){ // Electrons ...
3309  if(((positiveMC->GetMCProcessCode())) == 5 && ((negativeMC->GetMCProcessCode())) == 5){ // ... From Conversion ...
3310  if(gammaMC0->GetPdgCode() == 22){ // ... with Gamma Mother
3311  gamma0MotherLabel=gammaMC0->GetMother();
3312  }
3313  }
3314  if(gammaMC0->GetPdgCode() ==fMesonPDG){ // Dalitz candidate
3315  gamma0DalitzCand = kTRUE;
3316  gamma0MotherLabel=-fMesonPDG;
3317  }
3318  }
3319  }
3320  positiveMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate1->GetMCLabelPositive()));
3321  negativeMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate1->GetMCLabelNegative()));
3322 
3323  Int_t gamma1MCLabel = -1;
3324  Int_t gamma1MotherLabel = -1;
3325  if(!positiveMC||!negativeMC)
3326  return;
3327 
3328  if(positiveMC->GetMother()>-1&&(negativeMC->GetMother() == positiveMC->GetMother())){
3329  gamma1MCLabel = positiveMC->GetMother();
3330  }
3331  if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
3332  // Daughters Gamma 1
3333  AliAODMCParticle * gammaMC1 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MCLabel));
3334  if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){ // Electrons ...
3335  if(((positiveMC->GetMCProcessCode())) == 5 && ((negativeMC->GetMCProcessCode())) == 5){ // ... From Conversion ...
3336  if(gammaMC1->GetPdgCode() == 22){ // ... with Gamma Mother
3337  gamma1MotherLabel=gammaMC1->GetMother();
3338  }
3339  }
3340  if(gammaMC1->GetPdgCode() ==fMesonPDG ){ // Dalitz candidate
3341  gamma1DalitzCand = kTRUE;
3342  gamma1MotherLabel=-fMesonPDG;
3343  }
3344  }
3345  }
3346  if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
3347  if(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->GetPdgCode() == fMesonPDG){
3348  isTrueMeson=kTRUE;
3349  if (CheckVectorForDoubleCount(fVectorDoubleCountTrueMesons,gamma0MotherLabel)){
3350  if (fHistoDoubleCountTrueMesonInvMassPt[fiCut]) fHistoDoubleCountTrueMesonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),fWeightJetJetMC);
3351  FillMultipleCountMap(fMapMultipleCountTrueMesons,gamma0MotherLabel);
3352  }
3353  }
3354  }
3355 
3356  //Identify Dalitz candidate
3357  if (gamma1DalitzCand || gamma0DalitzCand){
3358  if (gamma0DalitzCand && gamma0MCLabel >=0 && gamma0MCLabel==gamma1MotherLabel){
3359  if (gamma0MotherLabel == -fMesonPDG) isTrueMesonDalitz = kTRUE;
3360  }
3361  if (gamma1DalitzCand && gamma1MCLabel >=0 && gamma1MCLabel==gamma0MotherLabel){
3362  if (gamma1MotherLabel == -fMesonPDG) isTrueMesonDalitz = kTRUE;
3363  }
3364  }
3365 
3366  if(isTrueMeson ){// True Meson
3367  if (fHistoTrueMesonInvMassPt[fiCut])fHistoTrueMesonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),fWeightJetJetMC);
3368  if (fDoMesonQA > 0){
3369  if ( mesonCand->M() > fMesonInvMassWindow[0] && mesonCand->M() < fMesonInvMassWindow[1]){
3370  if(fIsMC < 2){
3371  if (fHistoTrueMesonPtY[fiCut]) fHistoTrueMesonPtY[fiCut]->Fill(mesonCand->Pt(),mesonCand->Rapidity()-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift());
3372  if (fHistoTrueMesonPtOpenAngle[fiCut]) fHistoTrueMesonPtOpenAngle[fiCut]->Fill(mesonCand->Pt(),mesonCand->GetOpeningAngle());
3373  }
3374  if (fHistoTrueMesonPtAlpha[fiCut]) fHistoTrueMesonPtAlpha[fiCut]->Fill(mesonCand->Pt(),TMath::Abs(mesonCand->GetAlpha()),fWeightJetJetMC);
3375  }
3376 
3377  }
3378  Bool_t isPrimary = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsConversionPrimaryAOD(fInputEvent, static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma0MotherLabel)), mcProdVtxX, mcProdVtxY, mcProdVtxZ);
3379 
3380  if(isPrimary){ // Only primary pi0 for efficiency calculation
3381  Float_t weighted= 1;
3382  if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(gamma1MotherLabel, 0x0, fInputEvent)){
3383  if (static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt()>0.005){
3384  weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(gamma1MotherLabel, 0x0, fInputEvent);
3385  // cout << "rec \t " <<gamma1MotherLabel << "\t" << weighted << endl;
3386  }
3387  }
3388  if (fHistoTruePrimaryMesonInvMassPt[fiCut]) fHistoTruePrimaryMesonInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),weighted*fWeightJetJetMC);
3389  if (fHistoTruePrimaryMesonW0WeightingInvMassPt[fiCut]) fHistoTruePrimaryMesonW0WeightingInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),fWeightJetJetMC);
3390  if (fProfileTruePrimaryMesonWeightsInvMassPt[fiCut]) fProfileTruePrimaryMesonWeightsInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt(),weighted*fWeightJetJetMC);
3391 
3392  if (fDoMesonQA > 0 && fIsMC < 2){
3393  if (fHistoTruePrimaryMesonMCPtResolPt[fiCut])fHistoTruePrimaryMesonMCPtResolPt[fiCut]->Fill(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),
3394  (mesonCand->Pt()-static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt())/static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),weighted);
3395  }
3396  }
3397  } else if(!isTrueMeson ) { // Background
3398  if (fDoMesonQA > 1 && fIsMC < 2){
3399  if(!(gamma0MotherLabel>-1 && gamma1MotherLabel>-1)){ // Both Tracks are Photons and have a mother but not Pi0 or Eta
3400  if (fHistoTrueBckGGInvMassPt[fiCut])fHistoTrueBckGGInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
3401  } else { // No photon or without mother
3402  if (fHistoTrueBckContInvMassPt[fiCut]) fHistoTrueBckContInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
3403  }
3404  }
3405  if (!isTrueMesonDalitz && (gamma0DalitzCand || gamma1DalitzCand)){
3406  if (fDoMesonQA > 1 && fIsMC < 2)if (fHistoTrueBckContInvMassPt[fiCut]) fHistoTrueBckContInvMassPt[fiCut]->Fill(mesonCand->M(),mesonCand->Pt());
3407  }
3408  }
3409  }
3410  return;
3411 }
3412 
3413 
3414 //________________________________________________________________________
3415 void AliAnalysisTaskHeavyNeutralMesonToGG::CalculateBackground(){
3416 
3417  // set current BG handler
3418  AliGammaConversionAODBGHandler* currBGHandler = NULL;
3419  TList* currPhotonList = NULL;
3420  if (fMesonRecoMode == 0){
3421  currBGHandler = fBGHandler[fiCut];
3422  currPhotonList = fGammaCandidates;
3423  } else if (fMesonRecoMode == 1){
3424  currBGHandler = fBGClusHandler[fiCut];
3425  currPhotonList = fGammaCandidates;
3426  } else if (fMesonRecoMode == 2){
3427  currBGHandler = fBGClusHandler[fiCut];
3428  currPhotonList = fClusterCandidates;
3429  }
3430 
3431  Int_t zbin = currBGHandler->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
3432  Int_t mbin = 0;
3433  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
3434  mbin = currBGHandler->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
3435  } else {
3436  mbin = currBGHandler->GetMultiplicityBinIndex(currPhotonList->GetEntries());
3437  }
3438 
3439  AliGammaConversionAODBGHandler::GammaConversionVertex *bgEventVertex = NULL;
3440  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
3441  for(Int_t nEventsInBG=0;nEventsInBG<currBGHandler->GetNBGEvents();nEventsInBG++){
3442  AliGammaConversionAODVector *previousEventV0s = currBGHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);
3443  if(fMoveParticleAccordingToVertex == kTRUE || ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->GetInPlaneOutOfPlaneCut() != 0){
3444  bgEventVertex = currBGHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
3445  }
3446  for(Int_t iCurrent=0;iCurrent<currPhotonList->GetEntries();iCurrent++){
3447  AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(currPhotonList->At(iCurrent));
3448  for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
3449  AliAODConversionPhoton previousGoodV0 = (AliAODConversionPhoton)(*(previousEventV0s->at(iPrevious)));
3450  if(fMoveParticleAccordingToVertex == kTRUE){
3451  if (bgEventVertex){
3452  MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
3453  }
3454  }
3455  if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->GetInPlaneOutOfPlaneCut() != 0){
3456  if (bgEventVertex){
3457  RotateParticleAccordingToEP(&previousGoodV0,bgEventVertex->fEP,fEventPlaneAngle);
3458  }
3459  }
3460 
3461  AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&previousGoodV0);
3462  backgroundCandidate->CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
3463  if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift()))){
3464  if (fHistoMotherBackInvMassPt[fiCut]) fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt(),fWeightJetJetMC);
3465  if(fDoTHnSparse){
3466  Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin};
3467  if (fSparseMotherBackInvMassPtZM[fiCut]) fSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
3468  }
3469  }
3470  delete backgroundCandidate;
3471  backgroundCandidate = 0x0;
3472  }
3473  }
3474  }
3475  }else {
3476  for(Int_t nEventsInBG=0;nEventsInBG <currBGHandler->GetNBGEvents();nEventsInBG++){
3477  AliGammaConversionAODVector *previousEventV0s = currBGHandler->GetBGGoodV0s(zbin,mbin,nEventsInBG);
3478  if(previousEventV0s){
3479  if(fMoveParticleAccordingToVertex == kTRUE || ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->GetInPlaneOutOfPlaneCut() != 0 ){
3480  bgEventVertex = currBGHandler->GetBGEventVertex(zbin,mbin,nEventsInBG);
3481  }
3482  for(Int_t iCurrent=0;iCurrent<currPhotonList->GetEntries();iCurrent++){
3483  AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(currPhotonList->At(iCurrent));
3484  for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
3485 
3486  AliAODConversionPhoton previousGoodV0 = (AliAODConversionPhoton)(*(previousEventV0s->at(iPrevious)));
3487 
3488  if(fMoveParticleAccordingToVertex == kTRUE){
3489  if (bgEventVertex){
3490  MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
3491  }
3492  }
3493  if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->GetInPlaneOutOfPlaneCut() != 0){
3494  if (bgEventVertex){
3495  RotateParticleAccordingToEP(&previousGoodV0,bgEventVertex->fEP,fEventPlaneAngle);
3496  }
3497  }
3498 
3499  AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&previousGoodV0);
3500  backgroundCandidate->CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
3501  if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift()))){
3502  if (fHistoMotherBackInvMassPt[fiCut]) fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt(),fWeightJetJetMC);
3503  if(fDoTHnSparse){
3504  Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin};
3505  if (fSparseMotherBackInvMassPtZM[fiCut]) fSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
3506  }
3507  }
3508  delete backgroundCandidate;
3509  backgroundCandidate = 0x0;
3510  }
3511  }
3512  }
3513  }
3514  }
3515 }
3516 
3517 //________________________________________________________________________
3518 void AliAnalysisTaskHeavyNeutralMesonToGG::CalculateBackgroundRP(){
3519 
3520  // set current BG handler
3521  AliConversionAODBGHandlerRP* currBGHandler = NULL;
3522  TList* currPhotonList = NULL;
3523  if (fMesonRecoMode == 0){
3524  currBGHandler = fBGHandlerRP[fiCut];
3525  currPhotonList = fGammaCandidates;
3526  } else if (fMesonRecoMode == 1){
3527  currBGHandler = fBGClusHandlerRP[fiCut];
3528  currPhotonList = fGammaCandidates;
3529  } else if (fMesonRecoMode == 2){
3530  currBGHandler = fBGClusHandlerRP[fiCut];
3531  currPhotonList = fClusterCandidates;
3532  }
3533 
3534  Int_t zbin = currBGHandler->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
3535  Int_t mbin = 0;
3536  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
3537  mbin = currBGHandler->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
3538  } else {
3539  mbin = currBGHandler->GetMultiplicityBinIndex(currPhotonList->GetEntries());
3540  }
3541 
3542  //Rotation Method
3543  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseRotationMethod()){
3544  // Correct for the number of rotations
3545  // BG is for rotation the same, except for factor NRotations
3546  Double_t weight=1./Double_t(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->GetNumberOfBGEvents());
3547  for(Int_t firstGammaIndex=0;firstGammaIndex<currPhotonList->GetEntries();firstGammaIndex++){
3548  AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(currPhotonList->At(firstGammaIndex));
3549  if (gamma0==NULL) continue;
3550  for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<currPhotonList->GetEntries();secondGammaIndex++){
3551  AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(currPhotonList->At(secondGammaIndex));
3552  if (gamma1 == NULL) continue;
3553  if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelected(gamma1,fInputEvent))continue;
3554  for(Int_t nRandom=0;nRandom<((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->GetNumberOfBGEvents();nRandom++){
3555  RotateParticle(gamma1);
3556  AliAODConversionMother backgroundCandidate(gamma0,gamma1);
3557  backgroundCandidate.CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
3558  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(&backgroundCandidate, kFALSE, ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift())){
3559  if (fHistoMotherBackInvMassPt[fiCut]) fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate.M(),backgroundCandidate.Pt(),fWeightJetJetMC);
3560  if(fDoTHnSparse){
3561  Double_t sparesFill[4] = {backgroundCandidate.M(),backgroundCandidate.Pt(),(Double_t)zbin,(Double_t)mbin};
3562  if (fSparseMotherBackInvMassPtZM[fiCut]) fSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,weight);
3563  }
3564  }
3565  }
3566  }
3567  }
3568 
3569  }else {
3570  // Do Event Mixing
3571  for(Int_t nEventsInBG=0;nEventsInBG <currBGHandler->GetNBGEvents(currPhotonList,fInputEvent);nEventsInBG++){
3572  AliGammaConversionPhotonVector *previousEventGammas = currBGHandler->GetBGGoodGammas(currPhotonList,fInputEvent,nEventsInBG);
3573  if(previousEventGammas){
3574  // test weighted background
3575  Double_t weight=1.0;
3576  // Correct for the number of eventmixing:
3577  // N gammas -> (N-1) + (N-2) +(N-3) ...+ (N-(N-1)) using sum formula sum(i)=N*(N-1)/2 -> N*(N-1)/2
3578  // real combinations (since you cannot combine a photon with its own)
3579  // but BG leads to N_{a}*N_{b}combinations
3580  weight*=0.5*(Double_t(currPhotonList->GetEntries()-1))/Double_t(previousEventGammas->size());
3581 
3582  for(Int_t iCurrent=0;iCurrent<currPhotonList->GetEntries();iCurrent++){
3583  AliAODConversionPhoton *gamma0 = (AliAODConversionPhoton*)(currPhotonList->At(iCurrent));
3584  for(UInt_t iPrevious=0;iPrevious<previousEventGammas->size();iPrevious++){
3585  AliAODConversionPhoton *gamma1 = (AliAODConversionPhoton*)(previousEventGammas->at(iPrevious));
3586  AliAODConversionMother backgroundCandidate(gamma0,gamma1);
3587  backgroundCandidate.CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
3588  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
3589  ->MesonIsSelected(&backgroundCandidate,kFALSE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift())){
3590  if (fHistoMotherBackInvMassPt[fiCut]) fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate.M(),backgroundCandidate.Pt(),fWeightJetJetMC);
3591  if(fDoTHnSparse){
3592  Double_t sparesFill[4] = {backgroundCandidate.M(),backgroundCandidate.Pt(),(Double_t)zbin,(Double_t)mbin};
3593  if (fSparseMotherBackInvMassPtZM[fiCut]) fSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,weight);
3594  }
3595  }
3596  }
3597  }
3598  }
3599  }
3600  }
3601 }
3602 
3603 //________________________________________________________________________
3604 void AliAnalysisTaskHeavyNeutralMesonToGG::RotateParticle(AliAODConversionPhoton *gamma){
3605  Int_t fNDegreesPMBackground= ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->NDegreesRotation();
3606  Double_t nRadiansPM = fNDegreesPMBackground*TMath::Pi()/180;
3607  Double_t rotationValue = fRandom.Rndm()*2*nRadiansPM + TMath::Pi()-nRadiansPM;
3608  gamma->RotateZ(rotationValue);
3609 }
3610 
3611 //________________________________________________________________________
3612 void AliAnalysisTaskHeavyNeutralMesonToGG::RotateParticleAccordingToEP(AliAODConversionPhoton *gamma, Double_t previousEventEP, Double_t thisEventEP){
3613  previousEventEP=previousEventEP+TMath::Pi();
3614  thisEventEP=thisEventEP+TMath::Pi();
3615  Double_t rotationValue= thisEventEP-previousEventEP;
3616  gamma->RotateZ(rotationValue);
3617 }
3618 
3619 //________________________________________________________________________
3620 void AliAnalysisTaskHeavyNeutralMesonToGG::MoveParticleAccordingToVertex(AliAODConversionPhoton* particle,const AliGammaConversionAODBGHandler::GammaConversionVertex *vertex){
3621  //see header file for documentation
3622 
3623  Double_t dx = vertex->fX - fInputEvent->GetPrimaryVertex()->GetX();
3624  Double_t dy = vertex->fY - fInputEvent->GetPrimaryVertex()->GetY();
3625  Double_t dz = vertex->fZ - fInputEvent->GetPrimaryVertex()->GetZ();
3626 
3627  Double_t movedPlace[3] = {particle->GetConversionX() - dx,particle->GetConversionY() - dy,particle->GetConversionZ() - dz};
3628  particle->SetConversionPoint(movedPlace);
3629 }
3630 //________________________________________________________________________
3631 void AliAnalysisTaskHeavyNeutralMesonToGG::UpdateEventByEventData(){
3632  //see header file for documentation
3633  if(fGammaCandidates->GetEntries() >0 && fMesonRecoMode == 0 ){
3634  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
3635  fBGHandler[fiCut]->AddEvent(fGammaCandidates,fInputEvent->GetPrimaryVertex()->GetX(), fInputEvent->GetPrimaryVertex()->GetY(), fInputEvent->GetPrimaryVertex()->GetZ(),
3636  fV0Reader->GetNumberOfPrimaryTracks(), fEventPlaneAngle);
3637  }else { // means we use #V0s for multiplicity
3638  fBGHandler[fiCut]->AddEvent(fGammaCandidates, fInputEvent->GetPrimaryVertex()->GetX(), fInputEvent->GetPrimaryVertex()->GetY(), fInputEvent->GetPrimaryVertex()->GetZ(),
3639  fGammaCandidates->GetEntries(), fEventPlaneAngle);
3640  }
3641  } else if((fGammaCandidates->GetEntries() >0 || fClusterCandidates->GetEntries() > 0 ) && fMesonRecoMode == 1 ){
3642  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
3643  fBGHandler[fiCut]->AddEvent(fGammaCandidates,fInputEvent->GetPrimaryVertex()->GetX(), fInputEvent->GetPrimaryVertex()->GetY(), fInputEvent->GetPrimaryVertex()->GetZ(),
3644  fV0Reader->GetNumberOfPrimaryTracks(), fEventPlaneAngle);
3645  fBGClusHandler[fiCut]->AddEvent(fClusterCandidates,fInputEvent->GetPrimaryVertex()->GetX(), fInputEvent->GetPrimaryVertex()->GetY(), fInputEvent->GetPrimaryVertex()->GetZ(),
3646  fV0Reader->GetNumberOfPrimaryTracks(), fEventPlaneAngle);
3647  }else { // means we use #V0s for multiplicity
3648  fBGHandler[fiCut]->AddEvent(fGammaCandidates, fInputEvent->GetPrimaryVertex()->GetX(), fInputEvent->GetPrimaryVertex()->GetY(), fInputEvent->GetPrimaryVertex()->GetZ(),
3649  fGammaCandidates->GetEntries(), fEventPlaneAngle);
3650  fBGClusHandler[fiCut]->AddEvent(fClusterCandidates, fInputEvent->GetPrimaryVertex()->GetX(), fInputEvent->GetPrimaryVertex()->GetY(), fInputEvent->GetPrimaryVertex()->GetZ(),
3651  fGammaCandidates->GetEntries(), fEventPlaneAngle);
3652  }
3653  } else if(fClusterCandidates->GetEntries() >0 && fMesonRecoMode == 2 ){
3654  if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
3655  fBGClusHandler[fiCut]->AddEvent(fClusterCandidates,fInputEvent->GetPrimaryVertex()->GetX(), fInputEvent->GetPrimaryVertex()->GetY(), fInputEvent->GetPrimaryVertex()->GetZ(),
3656  fV0Reader->GetNumberOfPrimaryTracks(), fEventPlaneAngle);
3657  }else { // means we use #V0s for multiplicity
3658  fBGClusHandler[fiCut]->AddEvent(fClusterCandidates, fInputEvent->GetPrimaryVertex()->GetX(), fInputEvent->GetPrimaryVertex()->GetY(), fInputEvent->GetPrimaryVertex()->GetZ(),
3659  fGammaCandidates->GetEntries(), fEventPlaneAngle);
3660  }
3661  }
3662 }
3663 
3664 //________________________________________________________________________
3665 void AliAnalysisTaskHeavyNeutralMesonToGG::RelabelAODPhotonCandidates(Bool_t mode){
3666 
3667  // Relabeling For AOD Event
3668  // ESDiD -> AODiD
3669  // MCLabel -> AODMCLabel
3670 
3671  if(mode){
3672  fMCEventPos = new Int_t[fReaderGammas->GetEntries()];
3673  fMCEventNeg = new Int_t[fReaderGammas->GetEntries()];
3674  fESDArrayPos = new Int_t[fReaderGammas->GetEntries()];
3675  fESDArrayNeg = new Int_t[fReaderGammas->GetEntries()];
3676  }
3677 
3678  for(Int_t iGamma = 0;iGamma<fReaderGammas->GetEntries();iGamma++){
3679  AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) fReaderGammas->At(iGamma);
3680  if(!PhotonCandidate) continue;
3681  if(!mode){// Back to ESD Labels
3682  PhotonCandidate->SetMCLabelPositive(fMCEventPos[iGamma]);
3683  PhotonCandidate->SetMCLabelNegative(fMCEventNeg[iGamma]);
3684  PhotonCandidate->SetLabelPositive(fESDArrayPos[iGamma]);
3685  PhotonCandidate->SetLabelNegative(fESDArrayNeg[iGamma]);
3686  continue;
3687  }
3688  fMCEventPos[iGamma] = PhotonCandidate->GetMCLabelPositive();
3689  fMCEventNeg[iGamma] = PhotonCandidate->GetMCLabelNegative();
3690  fESDArrayPos[iGamma] = PhotonCandidate->GetTrackLabelPositive();
3691  fESDArrayNeg[iGamma] = PhotonCandidate->GetTrackLabelNegative();
3692 
3693  Bool_t AODLabelPos = kFALSE;
3694  Bool_t AODLabelNeg = kFALSE;
3695 
3696  for(Int_t i = 0; i<fInputEvent->GetNumberOfTracks();i++){
3697  AliAODTrack *tempDaughter = static_cast<AliAODTrack*>(fInputEvent->GetTrack(i));
3698  if(!AODLabelPos){
3699  if( tempDaughter->GetID() == PhotonCandidate->GetTrackLabelPositive() ){
3700  PhotonCandidate->SetMCLabelPositive(TMath::Abs(tempDaughter->GetLabel()));
3701  PhotonCandidate->SetLabelPositive(i);
3702  AODLabelPos = kTRUE;
3703  }
3704  }
3705  if(!AODLabelNeg){
3706  if( tempDaughter->GetID() == PhotonCandidate->GetTrackLabelNegative()){
3707  PhotonCandidate->SetMCLabelNegative(TMath::Abs(tempDaughter->GetLabel()));
3708  PhotonCandidate->SetLabelNegative(i);
3709  AODLabelNeg = kTRUE;
3710  }
3711  }
3712  if(AODLabelNeg && AODLabelPos){
3713  break;
3714  }
3715  }
3716  if(!AODLabelPos || !AODLabelNeg){
3717  cout<<"WARNING!!! AOD TRACKS NOT FOUND FOR"<<endl;
3718  if(!AODLabelNeg){
3719  PhotonCandidate->SetMCLabelNegative(-999999);
3720  PhotonCandidate->SetLabelNegative(-999999);
3721  }
3722  if(!AODLabelPos){
3723  PhotonCandidate->SetMCLabelPositive(-999999);
3724  PhotonCandidate->SetLabelPositive(-999999);
3725  }
3726  }
3727  }
3728 
3729  if(!mode){
3730  delete[] fMCEventPos;
3731  delete[] fMCEventNeg;
3732  delete[] fESDArrayPos;
3733  delete[] fESDArrayNeg;
3734  }
3735 }
3736 
3737 //________________________________________________________________________
3738 void AliAnalysisTaskHeavyNeutralMesonToGG::SetLogBinningXTH2(TH2* histoRebin){
3739  TAxis *axisafter = histoRebin->GetXaxis();
3740  Int_t bins = axisafter->GetNbins();
3741  Double_t from = axisafter->GetXmin();
3742  Double_t to = axisafter->GetXmax();
3743  Double_t *newbins = new Double_t[bins+1];
3744  newbins[0] = from;
3745  Double_t factor = TMath::Power(to/from, 1./bins);
3746  for(Int_t i=1; i<=bins; ++i) newbins[i] = factor * newbins[i-1];
3747  axisafter->Set(bins, newbins);
3748  delete [] newbins;
3749 }
3750 
3751 //________________________________________________________________________
3752 void AliAnalysisTaskHeavyNeutralMesonToGG::Terminate(const Option_t *)
3753 {
3754  //fOutputContainer->Print(); // Will crash on GRID
3755 }
3756 
3757 //_________________________________________________________________________________
3758 Bool_t AliAnalysisTaskHeavyNeutralMesonToGG::CheckVectorOnly(vector<Int_t> &vec, Int_t tobechecked){
3759  if(tobechecked > -1){
3760  vector<Int_t>::iterator it;
3761  it = find (vec.begin(), vec.end(), tobechecked);
3762  if (it != vec.end()) return true;
3763  else return false;
3764  }
3765  return false;
3766 }
3767 
3768 //_________________________________________________________________________________
3769 Bool_t AliAnalysisTaskHeavyNeutralMesonToGG::CheckVectorForDoubleCount(vector<Int_t> &vec, Int_t tobechecked){
3770  if(tobechecked > -1){
3771  vector<Int_t>::iterator it;
3772  it = find (vec.begin(), vec.end(), tobechecked);
3773  if (it != vec.end()) return true;
3774  else{
3775  vec.push_back(tobechecked);
3776  return false;
3777  }
3778  }
3779  return false;
3780 }
3781 
3782 //_________________________________________________________________________________
3783 void AliAnalysisTaskHeavyNeutralMesonToGG::FillMultipleCountMap(map<Int_t,Int_t> &ma, Int_t tobechecked){
3784  if( ma.find(tobechecked) != ma.end() ) ma[tobechecked] += 1;
3785  else ma[tobechecked] = 2;
3786  return;
3787 }
3788 
3789 //_________________________________________________________________________________
3790 void AliAnalysisTaskHeavyNeutralMesonToGG::FillMultipleCountHistoAndClear(map<Int_t,Int_t> &ma, TH1F* hist){
3791  map<Int_t, Int_t>::iterator it;
3792  for (it = ma.begin(); it != ma.end(); it++){
3793  hist->Fill(it->second, fWeightJetJetMC);
3794  }
3795  ma.clear();
3796  return;
3797 }
AliConversionPhotonBase::GetMCParticle
TParticle * GetMCParticle(AliMCEvent *mcEvent)
Definition: AliConversionPhotonBase.cxx:69
AliAnalysisTaskHeavyNeutralMesonToGG::fHistoClusOverlapHeadersGammaPt
TH1F ** fHistoClusOverlapHeadersGammaPt
array of histos with cluster, E
Definition: AliAnalysisTaskHeavyNeutralMesonToGG.h:197
AliAnalysisTaskHeavyNeutralMesonToGG::fHistoNGammaConvCandidates
TH1F ** fHistoNGammaConvCandidates
array of histos with vertex y distribution for selected events
Definition: AliAnalysisTaskHeavyNeutralMesonToGG.h:271
AliAODConversionPhoton::GetPhotonPhi
virtual Double_t GetPhotonPhi() const
Definition: AliAODConversionPhoton.h:41
AliAnalysisTaskHeavyNeutralMesonToGG::fHistoSPDClusterTrackletBackground
TH2F ** fHistoSPDClusterTrackletBackground
array of histos with number of calo gamma candidates per event
Definition: AliAnalysisTaskHeavyNeutralMesonToGG.h:273
AliAnalysisTaskEMCALClusterizeFast.h
AliAnalysisTaskHeavyNeutralMesonToGG::fHistoVertexX
TH1F ** fHistoVertexX
array of histos with vertex z distribution for selected events
Definition: AliAnalysisTaskHeavyNeutralMesonToGG.h:269
AliAnalysisTaskHeavyNeutralMesonToGG::fHistoTrueMesonCaloMergedClusterPartConvInvMassPt
TH2F ** fHistoTrueMesonCaloMergedClusterPartConvInvMassPt
array of histos with validated mothers, merged cluster invMass, pt
Definition: AliAnalysisTaskHeavyNeutralMesonToGG.h:224
AliAODConversionPhoton::SetCaloClusterRef
void SetCaloClusterRef(Long_t ref)
Definition: AliAODConversionPhoton.h:70
AliAODConversionMother::SetLabels
void SetLabels(Int_t label1, Int_t label2, Int_t label3=0)
Definition: AliAODConversionMother.h:58
Double_t
double Double_t
Definition: External.C:58
AliAnalysisTaskHeavyNeutralMesonToGG::CheckVectorForDoubleCount
Bool_t CheckVectorForDoubleCount(vector< Int_t > &vec, Int_t tobechecked)
Definition: AliAnalysisTaskHeavyNeutralMesonToGG.cxx:3769
AliAnalysisTaskHeavyNeutralMesonToGG.h
AliGammaConversionAODBGHandler::GetBGEventVertex
GammaConversionVertex * GetBGEventVertex(Int_t zbin, Int_t mbin, Int_t event)
Definition: AliGammaConversionAODBGHandler.h:81
AliAnalysisTaskHeavyNeutralMesonToGG::fHistoMCMesonWOWeightInAccPt
TH1F ** fHistoMCMesonWOWeightInAccPt
array of histos with weighted meson in acceptance, pT
Definition: AliAnalysisTaskHeavyNeutralMesonToGG.h:209
TH2F
Definition: External.C:236
AliAODConversionPhoton::SetCaloPhotonMCFlags
void SetCaloPhotonMCFlags(AliMCEvent *mcEvent, Bool_t enableSort)
Definition: AliAODConversionPhoton.cxx:147
AliConversionPhotonBase::SetConversionPoint
void SetConversionPoint(Double_t convpoint[3])
Definition: AliConversionPhotonBase.h:89