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