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