PWGGA/GammaConv/AliAnalysisTaskGammaConvCalo.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                                                                                                *
   4  * Author: Baldo Sahlmueller, Friederike Bock                                             *
   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
  17 //////////////////////////////////////////////////////////////////
  18 //----------------------------------------------------------------
  19 // Class used to do analysis on conversion photons + calo photons
  20 //----------------------------------------------------------------
  21 //////////////////////////////////////////////////////////////////
  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 "AliAnalysisTaskGammaConvCalo.h"
  43 #include "AliVParticle.h"
  44 #include "AliESDtrack.h"
  45 #include "AliESDtrackCuts.h"
  46 #include "AliKFVertex.h"
  47 #include "AliV0ReaderV1.h"
  48 #include "AliGenCocktailEventHeader.h"
  49 #include "AliConversionAODBGHandlerRP.h"
  50 #include "AliAODMCParticle.h"
  51 #include "AliAODMCHeader.h"
  52 #include "AliEventplane.h"
  53 #include "AliAnalysisTaskEMCALClusterizeFast.h"
  54 #include "AliAODEvent.h"
  55 #include "AliAODInputHandler.h"
  56 #include "AliESDEvent.h"
  57 #include "AliESDInputHandler.h"
  58 #include "AliInputEventHandler.h"
  59
  60 ClassImp(AliAnalysisTaskGammaConvCalo)
  61
  62 //________________________________________________________________________
  63 AliAnalysisTaskGammaConvCalo::AliAnalysisTaskGammaConvCalo(): AliAnalysisTaskSE(),
  64         fV0Reader(NULL),
  65         fBGHandler(NULL),
  66         fBGHandlerRP(NULL),
  67         fBGClusHandler(NULL),
  68         fBGClusHandlerRP(NULL),
  69         fInputEvent(NULL),
  70         fMCEvent(NULL),
  71         fMCStack(NULL),
  72         fCutFolder(NULL),
  73         fESDList(NULL),
  74         fBackList(NULL),
  75         fMotherList(NULL),
  76         fPhotonDCAList(NULL),
  77         fMesonDCAList(NULL),
  78         fTrueList(NULL),
  79         fMCList(NULL),
  80         fHeaderNameList(NULL),
  81         fTagOutputList(NULL),
  82         fOutputContainer(NULL),
  83         fReaderGammas(NULL),
  84         fGammaCandidates(NULL),
  85         fClusterCandidates(NULL),
  86         fEventCutArray(NULL),
  87         fEventCuts(NULL),
  88         fCutArray(NULL),
  89         fConversionCuts(NULL),
  90         fClusterCutArray(NULL),
  91         fCaloPhotonCuts(NULL),
  92         fMesonCutArray(NULL),
  93         fMesonCuts(NULL),
  94         fHistoConvGammaPt(NULL),
  95         fHistoConvGammaR(NULL),
  96         fHistoConvGammaEta(NULL),
  97         fTreeConvGammaPtDcazCat(NULL),
  98         fPtGamma(0),
  99         fDCAzPhoton(0),
 100         fRConvPhoton(0),
 101         fEtaPhoton(0),
 102         fCharCatPhoton(0),
 103         fCharPhotonMCInfo(0),
 104         fHistoMotherInvMassPt(NULL),
 105         fSparseMotherInvMassPtZM(NULL),
 106         fHistoMotherBackInvMassPt(NULL),
 107         fSparseMotherBackInvMassPtZM(NULL),
 108         fHistoMotherInvMassEalpha(NULL),
 109         fHistoMotherPi0PtY(NULL),
 110         fHistoMotherEtaPtY(NULL),
 111         fHistoMotherPi0PtAlpha(NULL),
 112         fHistoMotherEtaPtAlpha(NULL),
 113         fHistoMotherPi0PtOpenAngle(NULL),
 114         fHistoMotherEtaPtOpenAngle(NULL),
 115         fTreeMesonsInvMassPtDcazMinDcazMaxFlag(NULL),
 116         fInvMass(0),
 117         fPt(0),
 118         fDCAzGammaMin(0),
 119         fDCAzGammaMax(0),
 120         fCharFlag(0),
 121         fCharMesonMCInfo(0),
 122         fHistoConvGammaUntagged(NULL),
 123         fHistoConvGammaTagged(NULL),
 124         fHistoConvGammaPi0Tagged(NULL),
 125         fHistoConvGammaEtaTagged(NULL),
 126         fHistoPhotonPairAll(NULL),
 127         fHistoPhotonPairAllGam(NULL),
 128         fHistoClusGammaPt(NULL),
 129         fHistoMCHeaders(NULL),
 130         fHistoMCAllGammaPt(NULL),
 131         fHistoMCDecayGammaPi0Pt(NULL),
 132         fHistoMCDecayGammaRhoPt(NULL),
 133         fHistoMCDecayGammaEtaPt(NULL),
 134         fHistoMCDecayGammaOmegaPt(NULL),
 135         fHistoMCDecayGammaEtapPt(NULL),
 136         fHistoMCDecayGammaPhiPt(NULL),
 137         fHistoMCDecayGammaSigmaPt(NULL),
 138         fHistoMCConvGammaPt(NULL),
 139         fHistoMCConvGammaR(NULL),
 140         fHistoMCConvGammaEta(NULL),
 141         fHistoMCPi0Pt(NULL),
 142         fHistoMCPi0WOWeightPt(NULL),
 143         fHistoMCEtaPt(NULL),
 144         fHistoMCEtaWOWeightPt(NULL),
 145         fHistoMCPi0InAccPt(NULL),
 146         fHistoMCEtaInAccPt(NULL),
 147         fHistoMCPi0PtY(NULL),
 148         fHistoMCEtaPtY(NULL),
 149         fHistoMCK0sPt(NULL),
 150         fHistoMCK0sWOWeightPt(NULL),
 151         fHistoMCK0sPtY(NULL),
 152         fHistoMCSecPi0PtvsSource(NULL),
 153         fHistoMCSecPi0Source(NULL),
 154         fHistoMCSecEtaPt(NULL),
 155         fHistoMCSecEtaSource(NULL),
 156         fHistoTrueMotherInvMassPt(NULL),
 157         fHistoTruePrimaryMotherInvMassPt(NULL),
 158         fHistoTruePrimaryMotherW0WeightingInvMassPt(NULL),
 159         fProfileTruePrimaryMotherWeightsInvMassPt(NULL),
 160         fHistoTruePrimaryPi0MCPtResolPt(NULL),
 161         fHistoTruePrimaryEtaMCPtResolPt(NULL),
 162         fHistoTrueSecondaryMotherInvMassPt(NULL),
 163         fHistoTrueSecondaryMotherFromK0sInvMassPt(NULL),
 164         fHistoTrueK0sWithPi0DaughterMCPt(NULL),
 165         fHistoTrueSecondaryMotherFromEtaInvMassPt(NULL),
 166         fHistoTrueEtaWithPi0DaughterMCPt(NULL),
 167         fHistoTrueSecondaryMotherFromLambdaInvMassPt(NULL),
 168         fHistoTrueLambdaWithPi0DaughterMCPt(NULL),
 169         fHistoTrueBckGGInvMassPt(NULL),
 170         fHistoTrueBckContInvMassPt(NULL),
 171         fHistoTruePi0PtY(NULL),
 172         fHistoTrueEtaPtY(NULL),
 173         fHistoTruePi0PtAlpha(NULL),
 174         fHistoTrueEtaPtAlpha(NULL),
 175         fHistoTruePi0PtOpenAngle(NULL),
 176         fHistoTrueEtaPtOpenAngle(NULL),
 177         fHistoTrueMotherDalitzInvMassPt(NULL),
 178         fHistoTrueConvGammaPt(NULL),
 179         fHistoTrueConvPi0GammaPt(NULL),
 180         fHistoTrueConvGammaEta(NULL),
 181         fHistoCombinatorialPt(NULL),
 182         fHistoTruePrimaryConvGammaPt(NULL),
 183         fHistoTruePrimaryConvGammaESDPtMCPt(NULL),
 184         fHistoTrueSecondaryConvGammaPt(NULL),
 185         fHistoTrueSecondaryConvGammaFromXFromK0sPt(NULL),
 186         fHistoTrueSecondaryConvGammaFromXFromLambdaPt(NULL),
 187         fHistoTrueDalitzPsiPairDeltaPhi(NULL),
 188         fHistoTrueGammaPsiPairDeltaPhi(NULL),
 189         fHistoTrueClusGammaPt(NULL),
 190         fHistoTruePrimaryClusGammaPt(NULL),
 191         fHistoTruePrimaryClusGammaESDPtMCPt(NULL),
 192         fHistoNEvents(NULL),
 193         fHistoNGoodESDTracks(NULL),
 194         fHistoNGammaCandidates(NULL),
 195         fHistoNGoodESDTracksVsNGammaCanditates(NULL),
 196         fHistoNV0Tracks(NULL),
 197         fProfileEtaShift(NULL),
 198         fEventPlaneAngle(-100),
 199         fRandom(0),
 200         fNGammaCandidates(0),
 201         fUnsmearedPx(NULL),
 202         fUnsmearedPy(NULL),
 203         fUnsmearedPz(NULL),
 204         fUnsmearedE(NULL),
 205         fMCStackPos(NULL),
 206         fMCStackNeg(NULL),
 207         fESDArrayPos(NULL),
 208         fESDArrayNeg(NULL),
 209         fnCuts(0),
 210         fiCut(0),
 211         fMoveParticleAccordingToVertex(kTRUE),
 212         fIsHeavyIon(0),
 213         fDoMesonAnalysis(kTRUE),
 214         fDoMesonQA(0),
 215         fDoPhotonQA(0),
 216         fIsFromMBHeader(kTRUE),
 217         fIsMC(kFALSE),
 218         fMinE(0.1),
 219         fNminCells(2),
 220         fEMCm02cut(0.5)
 221 {
 222
 223 }
 224
 225 //________________________________________________________________________
 226 AliAnalysisTaskGammaConvCalo::AliAnalysisTaskGammaConvCalo(const char *name):
 227         AliAnalysisTaskSE(name),
 228         fV0Reader(NULL),
 229         fBGHandler(NULL),
 230         fBGHandlerRP(NULL),
 231         fBGClusHandler(NULL),
 232         fBGClusHandlerRP(NULL),
 233         fInputEvent(NULL),
 234         fMCEvent(NULL),
 235         fMCStack(NULL),
 236         fCutFolder(NULL),
 237         fESDList(NULL),
 238         fBackList(NULL),
 239         fMotherList(NULL),
 240         fPhotonDCAList(NULL),
 241         fMesonDCAList(NULL),
 242         fTrueList(NULL),
 243         fMCList(NULL),
 244         fHeaderNameList(NULL),
 245         fTagOutputList(NULL),
 246         fOutputContainer(0),
 247         fReaderGammas(NULL),
 248         fGammaCandidates(NULL),
 249         fClusterCandidates(NULL),
 250         fEventCutArray(NULL),
 251         fEventCuts(NULL),
 252         fCutArray(NULL),
 253         fConversionCuts(NULL),
 254         fClusterCutArray(NULL),
 255         fCaloPhotonCuts(NULL),
 256         fMesonCutArray(NULL),
 257         fMesonCuts(NULL),
 258         fHistoConvGammaPt(NULL),
 259         fHistoConvGammaR(NULL),
 260         fHistoConvGammaEta(NULL),
 261         fTreeConvGammaPtDcazCat(NULL),
 262         fPtGamma(0),
 263         fDCAzPhoton(0),
 264         fRConvPhoton(0),
 265         fEtaPhoton(0),
 266         fCharCatPhoton(0),
 267         fCharPhotonMCInfo(0),
 268         fHistoMotherInvMassPt(NULL),
 269         fSparseMotherInvMassPtZM(NULL),
 270         fHistoMotherBackInvMassPt(NULL),
 271         fSparseMotherBackInvMassPtZM(NULL),
 272         fHistoMotherInvMassEalpha(NULL),
 273         fHistoMotherPi0PtY(NULL),
 274         fHistoMotherEtaPtY(NULL),
 275         fHistoMotherPi0PtAlpha(NULL),
 276         fHistoMotherEtaPtAlpha(NULL),
 277         fHistoMotherPi0PtOpenAngle(NULL),
 278         fHistoMotherEtaPtOpenAngle(NULL),
 279         fTreeMesonsInvMassPtDcazMinDcazMaxFlag(NULL),
 280         fInvMass(0),
 281         fPt(0),
 282         fDCAzGammaMin(0),
 283         fDCAzGammaMax(0),
 284         fCharFlag(0),
 285         fCharMesonMCInfo(0),
 286         fHistoConvGammaUntagged(NULL),
 287         fHistoConvGammaTagged(NULL),
 288         fHistoConvGammaPi0Tagged(NULL),
 289         fHistoConvGammaEtaTagged(NULL),
 290         fHistoPhotonPairAll(NULL),
 291         fHistoPhotonPairAllGam(NULL),
 292         fHistoClusGammaPt(NULL),
 293         fHistoMCHeaders(NULL),
 294         fHistoMCAllGammaPt(NULL),
 295         fHistoMCDecayGammaPi0Pt(NULL),
 296         fHistoMCDecayGammaRhoPt(NULL),
 297         fHistoMCDecayGammaEtaPt(NULL),
 298         fHistoMCDecayGammaOmegaPt(NULL),
 299         fHistoMCDecayGammaEtapPt(NULL),
 300         fHistoMCDecayGammaPhiPt(NULL),
 301         fHistoMCDecayGammaSigmaPt(NULL),
 302         fHistoMCConvGammaPt(NULL),
 303         fHistoMCConvGammaR(NULL),
 304         fHistoMCConvGammaEta(NULL),
 305         fHistoMCPi0Pt(NULL),
 306         fHistoMCPi0WOWeightPt(NULL),
 307         fHistoMCEtaPt(NULL),
 308         fHistoMCEtaWOWeightPt(NULL),
 309         fHistoMCPi0InAccPt(NULL),
 310         fHistoMCEtaInAccPt(NULL),
 311         fHistoMCPi0PtY(NULL),
 312         fHistoMCEtaPtY(NULL),
 313         fHistoMCK0sPt(NULL),
 314         fHistoMCK0sWOWeightPt(NULL),
 315         fHistoMCK0sPtY(NULL),
 316         fHistoMCSecPi0PtvsSource(NULL),
 317         fHistoMCSecPi0Source(NULL),
 318         fHistoMCSecEtaPt(NULL),
 319         fHistoMCSecEtaSource(NULL),
 320         fHistoTrueMotherInvMassPt(NULL),
 321         fHistoTruePrimaryMotherInvMassPt(NULL),
 322         fHistoTruePrimaryMotherW0WeightingInvMassPt(NULL),
 323         fProfileTruePrimaryMotherWeightsInvMassPt(NULL),
 324         fHistoTruePrimaryPi0MCPtResolPt(NULL),
 325         fHistoTruePrimaryEtaMCPtResolPt(NULL),
 326         fHistoTrueSecondaryMotherInvMassPt(NULL),
 327         fHistoTrueSecondaryMotherFromK0sInvMassPt(NULL),
 328         fHistoTrueK0sWithPi0DaughterMCPt(NULL),
 329         fHistoTrueSecondaryMotherFromEtaInvMassPt(NULL),
 330         fHistoTrueEtaWithPi0DaughterMCPt(NULL),
 331         fHistoTrueSecondaryMotherFromLambdaInvMassPt(NULL),
 332         fHistoTrueLambdaWithPi0DaughterMCPt(NULL),
 333         fHistoTrueBckGGInvMassPt(NULL),
 334         fHistoTrueBckContInvMassPt(NULL),
 335         fHistoTruePi0PtY(NULL),
 336         fHistoTrueEtaPtY(NULL),
 337         fHistoTruePi0PtAlpha(NULL),
 338         fHistoTrueEtaPtAlpha(NULL),
 339         fHistoTruePi0PtOpenAngle(NULL),
 340         fHistoTrueEtaPtOpenAngle(NULL),
 341         fHistoTrueMotherDalitzInvMassPt(NULL),
 342         fHistoTrueConvGammaPt(NULL),
 343         fHistoTrueConvPi0GammaPt(NULL),
 344         fHistoTrueConvGammaEta(NULL),
 345         fHistoCombinatorialPt(NULL),
 346         fHistoTruePrimaryConvGammaPt(NULL),
 347         fHistoTruePrimaryConvGammaESDPtMCPt(NULL),
 348         fHistoTrueSecondaryConvGammaPt(NULL),
 349         fHistoTrueSecondaryConvGammaFromXFromK0sPt(NULL),
 350         fHistoTrueSecondaryConvGammaFromXFromLambdaPt(NULL),
 351         fHistoTrueDalitzPsiPairDeltaPhi(NULL),
 352         fHistoTrueGammaPsiPairDeltaPhi(NULL),
 353         fHistoTrueClusGammaPt(NULL),
 354         fHistoTruePrimaryClusGammaPt(NULL),
 355         fHistoTruePrimaryClusGammaESDPtMCPt(NULL),
 356         fHistoNEvents(NULL),
 357         fHistoNGoodESDTracks(NULL),
 358         fHistoNGammaCandidates(NULL),
 359         fHistoNGoodESDTracksVsNGammaCanditates(NULL),
 360         fHistoNV0Tracks(NULL),
 361         fProfileEtaShift(NULL),
 362         fEventPlaneAngle(-100),
 363         fRandom(0),
 364         fNGammaCandidates(0),
 365         fUnsmearedPx(NULL),
 366         fUnsmearedPy(NULL),
 367         fUnsmearedPz(NULL),
 368         fUnsmearedE(NULL),
 369         fMCStackPos(NULL),
 370         fMCStackNeg(NULL),
 371         fESDArrayPos(NULL),
 372         fESDArrayNeg(NULL),
 373         fnCuts(0),
 374         fiCut(0),
 375         fMoveParticleAccordingToVertex(kTRUE),
 376         fIsHeavyIon(0),
 377         fDoMesonAnalysis(kTRUE),
 378         fDoMesonQA(0),
 379         fDoPhotonQA(0),
 380         fIsFromMBHeader(kTRUE),
 381         fIsMC(kFALSE),
 382         fMinE(0.1),
 383         fNminCells(2),
 384         fEMCm02cut(0.5)
 385 {
 386   // Define output slots here
 387   DefineOutput(1, TList::Class());
 388 }
 389
 390 AliAnalysisTaskGammaConvCalo::~AliAnalysisTaskGammaConvCalo()
 391 {
 392         if(fGammaCandidates){
 393                 delete fGammaCandidates;
 394                 fGammaCandidates = 0x0;
 395         }
 396         if(fClusterCandidates){
 397                 delete fClusterCandidates;
 398                 fClusterCandidates = 0x0;
 399         }
 400         if(fBGHandler){
 401                 delete[] fBGHandler;
 402                 fBGHandler = 0x0;
 403         }
 404         if(fBGHandlerRP){
 405                 delete[] fBGHandlerRP;
 406                 fBGHandlerRP = 0x0;
 407         }
 408         if(fBGClusHandler){
 409                 delete[] fBGClusHandler;
 410                 fBGClusHandler = 0x0;
 411         }
 412         if(fBGClusHandlerRP){
 413                 delete[] fBGClusHandlerRP;
 414                 fBGClusHandlerRP = 0x0;
 415         }
 416 }
 417 //___________________________________________________________
 418 void AliAnalysisTaskGammaConvCalo::InitBack(){
 419
 420         const Int_t nDim = 4;
 421         Int_t nBins[nDim] = {800,250,7,4};
 422         Double_t xMin[nDim] = {0,0, 0,0};
 423         Double_t xMax[nDim] = {0.8,25,7,4};
 424
 425         fSparseMotherInvMassPtZM = new THnSparseF*[fnCuts];
 426         fSparseMotherBackInvMassPtZM = new THnSparseF*[fnCuts];
 427
 428         fBGHandler = new AliGammaConversionAODBGHandler*[fnCuts];
 429         fBGHandlerRP = new AliConversionAODBGHandlerRP*[fnCuts];
 430
 431         fBGClusHandler = new AliGammaConversionAODBGHandler*[fnCuts];
 432         fBGClusHandlerRP = new AliConversionAODBGHandlerRP*[fnCuts];
 433
 434         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 435                 if (((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->DoBGCalculation()){
 436                         TString cutstringEvent  = ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber();
 437                         TString cutstringPhoton = ((AliConversionPhotonCuts*)fCutArray->At(iCut))->GetCutNumber();
 438                         TString cutstringCalo   = ((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutNumber();
 439                         TString cutstringMeson  = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
 440
 441                         Int_t collisionSystem = atoi((TString)(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber())(0,1));
 442                         Int_t centMin = atoi((TString)(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber())(1,1));
 443                         Int_t centMax = atoi((TString)(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber())(2,1));
 444
 445                         if(collisionSystem == 1 || collisionSystem == 2 ||
 446                                 collisionSystem == 5 || collisionSystem == 8 ||
 447                                 collisionSystem == 9){
 448                                 centMin = centMin*10;
 449                                 centMax = centMax*10;
 450                                 if(centMax ==0 && centMax!=centMin) centMax=100;
 451                         } else if(collisionSystem == 3 || collisionSystem == 6){
 452                                 centMin = centMin*5;
 453                                 centMax = centMax*5;
 454                         } else if(collisionSystem == 4 || collisionSystem == 7){
 455                                 centMin = ((centMin*5)+45);
 456                                 centMax = ((centMax*5)+45);
 457                         }
 458
 459                         fBackList[iCut] = new TList();
 460                         fBackList[iCut]->SetName(Form("%s_%s_%s_%s Back histograms",cutstringEvent.Data(),cutstringPhoton.Data(),cutstringCalo.Data(), cutstringMeson.Data()));
 461                         fBackList[iCut]->SetOwner(kTRUE);
 462                         fCutFolder[iCut]->Add(fBackList[iCut]);
 463
 464                         fSparseMotherBackInvMassPtZM[iCut] = new THnSparseF("Back_Back_InvMass_Pt_z_m","Back_Back_InvMass_Pt_z_m",nDim,nBins,xMin,xMax);
 465                         fBackList[iCut]->Add(fSparseMotherBackInvMassPtZM[iCut]);
 466
 467                         fMotherList[iCut] = new TList();
 468                         fMotherList[iCut]->SetName(Form("%s_%s_%s_%s Mother histograms",cutstringEvent.Data(),cutstringPhoton.Data(),cutstringCalo.Data(),cutstringMeson.Data()));
 469                         fMotherList[iCut]->SetOwner(kTRUE);
 470                         fCutFolder[iCut]->Add(fMotherList[iCut]);
 471
 472                         fSparseMotherInvMassPtZM[iCut] = new THnSparseF("Back_Mother_InvMass_Pt_z_m","Back_Mother_InvMass_Pt_z_m",nDim,nBins,xMin,xMax);
 473                         fMotherList[iCut]->Add(fSparseMotherInvMassPtZM[iCut]);
 474
 475                         if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->BackgroundHandlerType() == 0){
 476                                 fBGHandler[iCut] = new AliGammaConversionAODBGHandler(
 477                                                                                                                                         collisionSystem,centMin,centMax,
 478                                                                                                                                         ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetNumberOfBGEvents(),
 479                                                                                                                                         ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseTrackMultiplicity());
 480                                 fBGClusHandler[iCut] = new AliGammaConversionAODBGHandler(
 481                                                                                                                                         collisionSystem,centMin,centMax,
 482                                                                                                                                         ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetNumberOfBGEvents(),
 483                                                                                                                                         ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseTrackMultiplicity());
 484                                 fBGHandlerRP[iCut] = NULL;
 485                         } else{
 486                                 fBGHandlerRP[iCut] = new AliConversionAODBGHandlerRP(
 487                                                                                                                                         ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsHeavyIon(),
 488                                                                                                                                         ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity(),
 489                                                                                                                                         ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetNumberOfBGEvents());
 490                                 fBGClusHandlerRP[iCut] = new AliConversionAODBGHandlerRP(
 491                                                                                                                                         ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsHeavyIon(),
 492                                                                                                                                         ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity(),
 493                                                                                                                                         ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetNumberOfBGEvents());
 494                                 fBGHandler[iCut] = NULL;
 495                         }
 496                 }
 497         }
 498 }
 499 //________________________________________________________________________
 500 void AliAnalysisTaskGammaConvCalo::UserCreateOutputObjects(){
 501
 502         // Create histograms
 503         if(fOutputContainer != NULL){
 504                 delete fOutputContainer;
 505                 fOutputContainer = NULL;
 506         }
 507         if(fOutputContainer == NULL){
 508                 fOutputContainer = new TList();
 509                 fOutputContainer->SetOwner(kTRUE);
 510         }
 511
 512         // Array of current cut's gammas
 513         fGammaCandidates = new TList();
 514         fClusterCandidates = new TList();
 515
 516         fCutFolder = new TList*[fnCuts];
 517         fESDList = new TList*[fnCuts];
 518         fBackList = new TList*[fnCuts];
 519         fMotherList = new TList*[fnCuts];
 520         fHistoNEvents = new TH1I*[fnCuts];
 521         fHistoNGoodESDTracks = new TH1I*[fnCuts];
 522         fHistoNGammaCandidates = new TH1I*[fnCuts];
 523         fHistoNGoodESDTracksVsNGammaCanditates = new TH2F*[fnCuts];
 524         fHistoNV0Tracks = new TH1I*[fnCuts];
 525         fProfileEtaShift = new TProfile*[fnCuts];
 526         fHistoConvGammaPt = new TH1F*[fnCuts];
 527
 528         if (fDoPhotonQA == 2){
 529                 fPhotonDCAList = new TList*[fnCuts];
 530                 fTreeConvGammaPtDcazCat = new TTree*[fnCuts];
 531         }
 532         if (fDoPhotonQA > 0){
 533                 fHistoConvGammaR = new TH1F*[fnCuts];
 534                 fHistoConvGammaEta = new TH1F*[fnCuts];
 535         }
 536
 537         if(fDoMesonAnalysis){
 538                 fHistoMotherInvMassPt = new TH2F*[fnCuts];
 539                 fHistoMotherBackInvMassPt = new TH2F*[fnCuts];
 540                 fHistoMotherInvMassEalpha = new TH2F*[fnCuts];
 541                 if (fDoMesonQA == 2){
 542                         fMesonDCAList = new TList*[fnCuts];
 543                         fTreeMesonsInvMassPtDcazMinDcazMaxFlag = new TTree*[fnCuts];
 544                 }
 545                 if (fDoMesonQA > 0){
 546                         fHistoMotherPi0PtY =  new TH2F*[fnCuts];
 547                         fHistoMotherEtaPtY =  new TH2F*[fnCuts];
 548                         fHistoMotherPi0PtAlpha =  new TH2F*[fnCuts];
 549                         fHistoMotherEtaPtAlpha =  new TH2F*[fnCuts];
 550                         fHistoMotherPi0PtOpenAngle =  new TH2F*[fnCuts];
 551                         fHistoMotherEtaPtOpenAngle =  new TH2F*[fnCuts];
 552                 }
 553         }
 554         fTagOutputList = new TList*[fnCuts];
 555
 556         fHistoConvGammaUntagged = new TH1F*[fnCuts];
 557         fHistoConvGammaTagged = new TH1F*[fnCuts];
 558         fHistoConvGammaPi0Tagged = new TH1F*[fnCuts];
 559         fHistoConvGammaEtaTagged = new TH1F*[fnCuts];
 560         fHistoPhotonPairAll = new TH2F*[fnCuts];
 561         fHistoPhotonPairAllGam = new TH2F*[fnCuts];
 562
 563         fHistoClusGammaPt = new TH1F*[fnCuts];
 564
 565         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 566                 TString cutstringEvent  = ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber();
 567                 TString cutstringPhoton = ((AliConversionPhotonCuts*)fCutArray->At(iCut))->GetCutNumber();
 568                 TString cutstringCalo   = ((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutNumber();
 569                 TString cutstringMeson  = "NoMesonCut";
 570                 if(fDoMesonAnalysis)cutstringMeson = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
 571
 572                 fCutFolder[iCut] = new TList();
 573                 fCutFolder[iCut]->SetName(Form("Cut Number %s_%s_%s_%s",cutstringEvent.Data(),cutstringPhoton.Data(),cutstringCalo.Data(),cutstringMeson.Data()));
 574                 fCutFolder[iCut]->SetOwner(kTRUE);
 575                 fOutputContainer->Add(fCutFolder[iCut]);
 576                 fESDList[iCut] = new TList();
 577                 fESDList[iCut]->SetName(Form("%s_%s_%s_%s ESD histograms",cutstringEvent.Data(),cutstringPhoton.Data(),cutstringCalo.Data(),cutstringMeson.Data()));
 578                 fESDList[iCut]->SetOwner(kTRUE);
 579                 fCutFolder[iCut]->Add(fESDList[iCut]);
 580
 581                 fHistoNEvents[iCut] = new TH1I("NEvents","NEvents",9,-0.5,8.5);
 582                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(1,"Accepted");
 583                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(2,"Centrality");
 584                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(3,"Missing MC");
 585                 if (((AliConvEventCuts*)fEventCutArray->At(iCut))->IsSpecialTrigger() == 4 ){
 586                         TString TriggerNames = "Not Trigger: ";
 587                         TriggerNames = TriggerNames+ ( (AliConvEventCuts*)fEventCutArray->At(iCut))->GetSpecialTriggerName();
 588                         fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(4,TriggerNames.Data());
 589                 } else {
 590                         fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(4,"Trigger");
 591                 }
 592                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(5,"Vertex Z");
 593                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(6,"Cont. Vertex");
 594                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(7,"Pile-Up");
 595                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(8,"no SDD");
 596                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(9,"no V0AND");
 597                 fESDList[iCut]->Add(fHistoNEvents[iCut]);
 598
 599                 if(fIsHeavyIon == 1) fHistoNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",4000,0,4000);
 600                 else if(fIsHeavyIon == 2) fHistoNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",400,0,400);
 601                 else fHistoNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",200,0,200);
 602                 fESDList[iCut]->Add(fHistoNGoodESDTracks[iCut]);
 603                 if(fIsHeavyIon == 1) fHistoNGammaCandidates[iCut] = new TH1I("GammaCandidates","GammaCandidates",100,0,100);
 604                 else if(fIsHeavyIon == 2) fHistoNGammaCandidates[iCut] = new TH1I("GammaCandidates","GammaCandidates",50,0,50);
 605                 else fHistoNGammaCandidates[iCut] = new TH1I("GammaCandidates","GammaCandidates",50,0,50);
 606                 fESDList[iCut]->Add(fHistoNGammaCandidates[iCut]);
 607                 if(fIsHeavyIon == 1) fHistoNGoodESDTracksVsNGammaCanditates[iCut] = new TH2F("GoodESDTracksVsGammaCandidates","GoodESDTracksVsGammaCandidates",4000,0,4000,100,0,100);
 608                 else if(fIsHeavyIon == 2) fHistoNGoodESDTracksVsNGammaCanditates[iCut] = new TH2F("GoodESDTracksVsGammaCandidates","GoodESDTracksVsGammaCandidates",400,0,400,50,0,50);
 609                 else fHistoNGoodESDTracksVsNGammaCanditates[iCut] = new TH2F("GoodESDTracksVsGammaCandidates","GoodESDTracksVsGammaCandidates",200,0,200,50,0,50);
 610                 fESDList[iCut]->Add(fHistoNGoodESDTracksVsNGammaCanditates[iCut]);
 611
 612
 613                 if(fIsHeavyIon == 1) fHistoNV0Tracks[iCut] = new TH1I("V0 Multiplicity","V0 Multiplicity",30000,0,30000);
 614                 else if(fIsHeavyIon == 2) fHistoNV0Tracks[iCut] = new TH1I("V0 Multiplicity","V0 Multiplicity",2500,0,2500);
 615                 else fHistoNV0Tracks[iCut] = new TH1I("V0 Multiplicity","V0 Multiplicity",1500,0,1500);
 616                 fESDList[iCut]->Add(fHistoNV0Tracks[iCut]);
 617                 fProfileEtaShift[iCut] = new TProfile("Eta Shift","Eta Shift",1, -0.5,0.5);
 618                 fESDList[iCut]->Add(fProfileEtaShift[iCut]);
 619                 fHistoConvGammaPt[iCut] = new TH1F("ESD_ConvGamma_Pt","ESD_ConvGamma_Pt",250,0,25);
 620                 fESDList[iCut]->Add(fHistoConvGammaPt[iCut]);
 621
 622                 if (fDoPhotonQA == 2){
 623                         fPhotonDCAList[iCut] = new TList();
 624                         fPhotonDCAList[iCut]->SetName(Form("%s_%s_%s_%s Photon DCA tree",cutstringEvent.Data(),cutstringPhoton.Data(),cutstringCalo.Data(),cutstringMeson.Data()));
 625                         fPhotonDCAList[iCut]->SetOwner(kTRUE);
 626                         fCutFolder[iCut]->Add(fPhotonDCAList[iCut]);
 627
 628                         fTreeConvGammaPtDcazCat[iCut] = new TTree("ESD_ConvGamma_Pt_Dcaz_R_Eta","ESD_ConvGamma_Pt_Dcaz_R_Eta_Cat");
 629                         fTreeConvGammaPtDcazCat[iCut]->Branch("Pt",&fPtGamma,"fPtGamma/F");
 630                         fTreeConvGammaPtDcazCat[iCut]->Branch("DcaZPhoton",&fDCAzPhoton,"fDCAzPhoton/F");
 631       //          fTreeConvGammaPtDcazCat[iCut]->Branch("R",&fRConvPhoton,"fRConvPhoton/F");
 632       //          fTreeConvGammaPtDcazCat[iCut]->Branch("Eta",&fEtaPhoton,"fEtaPhoton/F");
 633
 634                         fTreeConvGammaPtDcazCat[iCut]->Branch("cat",&fCharCatPhoton,"fCharCatPhoton/b");
 635                         if(fIsMC){
 636                                 fTreeConvGammaPtDcazCat[iCut]->Branch("photonMCInfo",&fCharPhotonMCInfo,"fCharPhotonMCInfo/b");
 637                         }
 638                         fPhotonDCAList[iCut]->Add(fTreeConvGammaPtDcazCat[iCut]);
 639                 }
 640
 641                 if (fDoPhotonQA > 0){
 642                         fHistoConvGammaR[iCut] = new TH1F("ESD_ConvGamma_R","ESD_ConvGamma_R",800,0,200);
 643                         fESDList[iCut]->Add(fHistoConvGammaR[iCut]);
 644                         fHistoConvGammaEta[iCut] = new TH1F("ESD_ConvGamma_Eta","ESD_ConvGamma_Eta",2000,-2,2);
 645                         fESDList[iCut]->Add(fHistoConvGammaEta[iCut]);
 646                 }
 647
 648                 fTagOutputList[iCut] = new TList();
 649                 fTagOutputList[iCut]->SetName(Form("%s_%s_%s_%s Tagging Output",cutstringEvent.Data(),cutstringPhoton.Data(),cutstringCalo.Data(),cutstringMeson.Data()));
 650                 fTagOutputList[iCut]->SetOwner(1);
 651                 fCutFolder[iCut]->Add(fTagOutputList[iCut]);
 652
 653                 const Int_t nptbins = 200;
 654                 const Double_t ptmin = 0.;
 655                 const Double_t ptmax = 20.;
 656
 657                 const Int_t nmbins = 180;
 658                 const Double_t mmin = 0.;
 659                 const Double_t mmax = 0.9;
 660
 661                 // photon candidates
 662                 // this is maybe not necessary ...
 663
 664                 fHistoConvGammaUntagged[iCut] = new TH1F("ConvGammaUntagged","",nptbins,ptmin,ptmax);
 665                 fHistoConvGammaUntagged[iCut]->SetXTitle("p_{T} (GeV/c)");
 666                 fTagOutputList[iCut]->Add(fHistoConvGammaUntagged[iCut]);
 667
 668                 fHistoConvGammaTagged[iCut] = new TH1F("ConvGammaTagged","",nptbins,ptmin,ptmax);
 669                 fHistoConvGammaTagged[iCut]->SetXTitle("p_{T} (GeV/c)");
 670                 fTagOutputList[iCut]->Add(fHistoConvGammaTagged[iCut]);
 671
 672                 fHistoConvGammaPi0Tagged[iCut] = new TH1F("ConvGammaPi0Tagged","",nptbins,ptmin,ptmax);
 673                 fHistoConvGammaPi0Tagged[iCut]->SetXTitle("p_{T} (GeV/c)");
 674                 fTagOutputList[iCut]->Add(fHistoConvGammaPi0Tagged[iCut]);
 675
 676                 fHistoConvGammaEtaTagged[iCut] = new TH1F("ConvGammaEtaTagged","",nptbins,ptmin,ptmax);
 677                 fHistoConvGammaEtaTagged[iCut]->SetXTitle("p_{T} (GeV/c)");
 678                 fTagOutputList[iCut]->Add(fHistoConvGammaEtaTagged[iCut]);
 679
 680                 // pairs
 681                 fHistoPhotonPairAll[iCut] = new TH2F("PhotonPairAll","",nmbins,mmin,mmax,nptbins,ptmin,ptmax);
 682                 fHistoPhotonPairAll[iCut]->SetXTitle("M_{inv} (GeV/cc)");
 683                 fHistoPhotonPairAll[iCut]->SetYTitle("p_{T} (GeV/c)");
 684                 fTagOutputList[iCut]->Add(fHistoPhotonPairAll[iCut]);
 685
 686                 fHistoPhotonPairAllGam[iCut] = new TH2F("PhotonPairAllGammaConvPt","",nmbins,mmin,mmax,nptbins,ptmin,ptmax);
 687                 fHistoPhotonPairAllGam[iCut]->SetXTitle("M_{inv} (GeV/cc)");
 688                 fHistoPhotonPairAllGam[iCut]->SetYTitle("#gamma^{conv} p_{T} (GeV/c)");
 689                 fTagOutputList[iCut]->Add(fHistoPhotonPairAllGam[iCut]);
 690
 691                 fHistoClusGammaPt[iCut] = new TH1F("ClusGamma_Pt","ClusGamma_Pt",250,0,25);
 692                 fTagOutputList[iCut]->Add(fHistoClusGammaPt[iCut]);
 693
 694
 695                 if(fDoMesonAnalysis){
 696                         fHistoMotherInvMassPt[iCut] = new TH2F("ESD_Mother_InvMass_Pt","ESD_Mother_InvMass_Pt",800,0,0.8,250,0,25);
 697                         fESDList[iCut]->Add(fHistoMotherInvMassPt[iCut]);
 698                         fHistoMotherBackInvMassPt[iCut] = new TH2F("ESD_Background_InvMass_Pt","ESD_Background_InvMass_Pt",800,0,0.8,250,0,25);
 699                         fESDList[iCut]->Add(fHistoMotherBackInvMassPt[iCut]);
 700                         fHistoMotherInvMassEalpha[iCut] = new TH2F("ESD_Mother_InvMass_vs_E_alpha","ESD_Mother_InvMass_vs_E_alpha",800,0,0.8,250,0,25);
 701                         fESDList[iCut]->Add(fHistoMotherInvMassEalpha[iCut]);
 702                         if (fDoMesonQA == 2){
 703                                 fMesonDCAList[iCut] = new TList();
 704                                 fMesonDCAList[iCut]->SetName(Form("%s_%s_%s_%s Meson DCA tree",cutstringEvent.Data(),cutstringPhoton.Data(),cutstringCalo.Data(),cutstringMeson.Data()));
 705                                 fMesonDCAList[iCut]->SetOwner(kTRUE);
 706                                 fCutFolder[iCut]->Add(fMesonDCAList[iCut]);
 707
 708                                 fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut] = new TTree("ESD_Mesons_InvMass_Pt_DcazMin_DcazMax_Flag","ESD_Mesons_InvMass_Pt_DcazMin_DcazMax_Flag");
 709                                 fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("InvMass",&fInvMass,"fInvMass/F");
 710                                 fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("Pt",&fPt,"fPt/F");
 711                                 fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("DcaZMin",&fDCAzGammaMin,"fDCAzGammaMin/F");
 712                                 fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("DcaZMax",&fDCAzGammaMax,"fDCAzGammaMax/F");
 713                                 fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("kind",&fCharFlag,"fCharFlag/b");
 714                                 if(fIsMC){
 715                                         fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("mesonMCInfo",&fCharMesonMCInfo,"fCharMesonMCInfo/b");
 716                                 }
 717                                 fMesonDCAList[iCut]->Add(fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]);
 718
 719                         }
 720                         if (fDoMesonQA > 0 ){
 721                                 fHistoMotherPi0PtY[iCut] = new TH2F("ESD_MotherPi0_Pt_Y","ESD_MotherPi0_Pt_Y",150,0.03,15.,150,-1.5,1.5);
 722                                 SetLogBinningXTH2(fHistoMotherPi0PtY[iCut]);
 723                                 fESDList[iCut]->Add(fHistoMotherPi0PtY[iCut]);
 724                                 fHistoMotherEtaPtY[iCut] = new TH2F("ESD_MotherEta_Pt_Y","ESD_MotherEta_Pt_Y",150,0.03,15.,150,-1.5,1.5);
 725                                 SetLogBinningXTH2(fHistoMotherEtaPtY[iCut]);
 726                                 fESDList[iCut]->Add(fHistoMotherEtaPtY[iCut]);
 727                                 fHistoMotherPi0PtAlpha[iCut] = new TH2F("ESD_MotherPi0_Pt_Alpha","ESD_MotherPi0_Pt_Alpha",150,0.03,15.,100,0,1);
 728                                 SetLogBinningXTH2(fHistoMotherPi0PtAlpha[iCut]);
 729                                 fESDList[iCut]->Add(fHistoMotherPi0PtAlpha[iCut]);
 730                                 fHistoMotherEtaPtAlpha[iCut] = new TH2F("ESD_MotherEta_Pt_Alpha","ESD_MotherEta_Pt_Alpha",150,0.03,15.,100,0,1);
 731                                 SetLogBinningXTH2(fHistoMotherEtaPtAlpha[iCut]);
 732                                 fESDList[iCut]->Add(fHistoMotherEtaPtAlpha[iCut]);
 733                                 fHistoMotherPi0PtOpenAngle[iCut] = new TH2F("ESD_MotherPi0_Pt_OpenAngle","ESD_MotherPi0_Pt_OpenAngle",150,0.03,15.,200,0,2*TMath::Pi());
 734                                 SetLogBinningXTH2(fHistoMotherPi0PtOpenAngle[iCut]);
 735                                 fESDList[iCut]->Add(fHistoMotherPi0PtOpenAngle[iCut]);
 736                                 fHistoMotherEtaPtOpenAngle[iCut] = new TH2F("ESD_MotherEta_Pt_OpenAngle","ESD_MotherEta_Pt_OpenAngle",150,0.03,15.,200,0,2*TMath::Pi());
 737                                 SetLogBinningXTH2(fHistoMotherEtaPtOpenAngle[iCut]);
 738                                 fESDList[iCut]->Add(fHistoMotherEtaPtOpenAngle[iCut]);
 739                         }
 740                 }
 741         }
 742         if(fDoMesonAnalysis){
 743                 InitBack(); // Init Background Handler
 744         }
 745
 746         if(fIsMC){
 747                 // MC Histogramms
 748                 fMCList = new TList*[fnCuts];
 749                 // True Histogramms
 750                 fTrueList = new TList*[fnCuts];
 751                 // Selected Header List
 752                 fHeaderNameList = new TList*[fnCuts];
 753                 fHistoMCHeaders = new TH1I*[fnCuts];
 754                 fHistoMCAllGammaPt = new TH1F*[fnCuts];
 755                 fHistoMCDecayGammaPi0Pt = new TH1F*[fnCuts];
 756                 fHistoMCDecayGammaRhoPt = new TH1F*[fnCuts];
 757                 fHistoMCDecayGammaEtaPt = new TH1F*[fnCuts];
 758                 fHistoMCDecayGammaOmegaPt = new TH1F*[fnCuts];
 759                 fHistoMCDecayGammaEtapPt = new TH1F*[fnCuts];
 760                 fHistoMCDecayGammaPhiPt = new TH1F*[fnCuts];
 761                 fHistoMCDecayGammaSigmaPt = new TH1F*[fnCuts];
 762                 fHistoMCConvGammaPt = new TH1F*[fnCuts];
 763                 fHistoTrueConvGammaPt = new TH1F*[fnCuts];
 764                 fHistoTrueConvPi0GammaPt = new TH1F*[fnCuts];
 765
 766                 fHistoCombinatorialPt = new TH2F*[fnCuts];
 767                 fHistoTruePrimaryConvGammaPt = new TH1F*[fnCuts];
 768                 fHistoTruePrimaryConvGammaESDPtMCPt = new TH2F*[fnCuts];
 769                 fHistoTrueSecondaryConvGammaPt = new TH1F*[fnCuts];
 770                 fHistoTrueSecondaryConvGammaFromXFromK0sPt = new TH1F*[fnCuts];
 771                 fHistoTrueSecondaryConvGammaFromXFromLambdaPt = new TH1F*[fnCuts];
 772
 773                 fHistoTrueDalitzPsiPairDeltaPhi= new TH2F*[fnCuts];
 774                 fHistoTrueGammaPsiPairDeltaPhi= new TH2F*[fnCuts];
 775
 776                 fHistoTrueClusGammaPt = new TH1F*[fnCuts];
 777                 fHistoTruePrimaryClusGammaPt = new TH1F*[fnCuts];
 778                 fHistoTruePrimaryClusGammaESDPtMCPt = new TH2F*[fnCuts];
 779
 780                 if (fDoPhotonQA > 0){
 781                         fHistoMCConvGammaR = new TH1F*[fnCuts];
 782                         fHistoMCConvGammaEta = new TH1F*[fnCuts];
 783                         fHistoTrueConvGammaEta = new TH1F*[fnCuts];
 784                 }
 785
 786                 if(fDoMesonAnalysis){
 787                         fHistoMCPi0Pt = new TH1F*[fnCuts];
 788                         fHistoMCPi0WOWeightPt = new TH1F*[fnCuts];
 789                         fHistoMCEtaPt = new TH1F*[fnCuts];
 790                         fHistoMCEtaWOWeightPt = new TH1F*[fnCuts];
 791                         fHistoMCPi0InAccPt = new TH1F*[fnCuts];
 792                         fHistoMCEtaInAccPt = new TH1F*[fnCuts];
 793
 794                         fHistoTrueMotherInvMassPt = new TH2F*[fnCuts];
 795                         fHistoTruePrimaryMotherInvMassPt = new TH2F*[fnCuts];
 796                         fHistoTruePrimaryMotherW0WeightingInvMassPt = new TH2F*[fnCuts];
 797                         fProfileTruePrimaryMotherWeightsInvMassPt = new TProfile2D*[fnCuts];
 798                         fHistoTrueSecondaryMotherInvMassPt = new TH2F*[fnCuts];
 799                         fHistoTrueSecondaryMotherFromK0sInvMassPt = new TH2F*[fnCuts];
 800                         fHistoTrueSecondaryMotherFromEtaInvMassPt = new TH2F*[fnCuts];
 801                         fHistoTrueSecondaryMotherFromLambdaInvMassPt = new TH2F*[fnCuts];
 802                         fHistoTrueMotherDalitzInvMassPt = new TH2F*[fnCuts];
 803                         if (fDoMesonQA > 0){
 804                                 fHistoMCPi0PtY = new TH2F*[fnCuts];
 805                                 fHistoMCEtaPtY = new TH2F*[fnCuts];
 806                                 fHistoMCK0sPt = new TH1F*[fnCuts];
 807                                 fHistoMCK0sWOWeightPt = new TH1F*[fnCuts];
 808                                 fHistoMCK0sPtY = new TH2F*[fnCuts];
 809                                 fHistoMCSecPi0PtvsSource= new TH2F*[fnCuts];
 810                                 fHistoMCSecPi0Source = new TH1F*[fnCuts];
 811                                 fHistoMCSecEtaPt = new TH1F*[fnCuts];
 812                                 fHistoMCSecEtaSource = new TH1F*[fnCuts];
 813                                 fHistoTruePrimaryPi0MCPtResolPt = new TH2F*[fnCuts];
 814                                 fHistoTruePrimaryEtaMCPtResolPt = new TH2F*[fnCuts];
 815                                 fHistoTrueK0sWithPi0DaughterMCPt = new TH1F*[fnCuts];
 816                                 fHistoTrueEtaWithPi0DaughterMCPt = new TH1F*[fnCuts];
 817                                 fHistoTrueLambdaWithPi0DaughterMCPt = new TH1F*[fnCuts];
 818                                 fHistoTrueBckGGInvMassPt = new TH2F*[fnCuts];
 819                                 fHistoTrueBckContInvMassPt = new TH2F*[fnCuts];
 820                                 fHistoTruePi0PtY = new TH2F*[fnCuts];
 821                                 fHistoTrueEtaPtY = new TH2F*[fnCuts];
 822                                 fHistoTruePi0PtAlpha = new TH2F*[fnCuts];
 823                                 fHistoTrueEtaPtAlpha = new TH2F*[fnCuts];
 824                                 fHistoTruePi0PtOpenAngle = new TH2F*[fnCuts];
 825                                 fHistoTrueEtaPtOpenAngle = new TH2F*[fnCuts];
 826                         }
 827                 }
 828
 829
 830
 831                 for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 832                         TString cutstringEvent  = ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber();
 833                         TString cutstringPhoton = ((AliConversionPhotonCuts*)fCutArray->At(iCut))->GetCutNumber();
 834                         TString cutstringCalo   = ((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutNumber();
 835                         TString cutstringMeson  = "NoMesonCut";
 836                         if(fDoMesonAnalysis)cutstringMeson = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
 837
 838                         fMCList[iCut] = new TList();
 839                         fMCList[iCut]->SetName(Form("%s_%s_%s_%s MC histograms",cutstringEvent.Data(),cutstringPhoton.Data(),cutstringCalo.Data(),cutstringMeson.Data()));
 840                         fMCList[iCut]->SetOwner(kTRUE);
 841                         fCutFolder[iCut]->Add(fMCList[iCut]);
 842                         fHistoMCHeaders[iCut] = new TH1I("MC_Headers","MC_Headers",20,0,20);
 843                         fMCList[iCut]->Add(fHistoMCHeaders[iCut]);
 844                         fHistoMCAllGammaPt[iCut] = new TH1F("MC_AllGamma_Pt","MC_AllGamma_Pt",250,0,25);
 845                         fMCList[iCut]->Add(fHistoMCAllGammaPt[iCut]);
 846                         fHistoMCDecayGammaPi0Pt[iCut] = new TH1F("MC_DecayGammaPi0_Pt","MC_DecayGammaPi0_Pt",250,0,25);
 847                         fMCList[iCut]->Add(fHistoMCDecayGammaPi0Pt[iCut]);
 848                         fHistoMCDecayGammaRhoPt[iCut] = new TH1F("MC_DecayGammaRho_Pt","MC_DecayGammaRho_Pt",250,0,25);
 849                         fMCList[iCut]->Add(fHistoMCDecayGammaRhoPt[iCut]);
 850                         fHistoMCDecayGammaEtaPt[iCut] = new TH1F("MC_DecayGammaEta_Pt","MC_DecayGammaEta_Pt",250,0,25);
 851                         fMCList[iCut]->Add(fHistoMCDecayGammaEtaPt[iCut]);
 852                         fHistoMCDecayGammaOmegaPt[iCut] = new TH1F("MC_DecayGammaOmega_Pt","MC_DecayGammaOmmega_Pt",250,0,25);
 853                         fMCList[iCut]->Add(fHistoMCDecayGammaOmegaPt[iCut]);
 854                         fHistoMCDecayGammaEtapPt[iCut] = new TH1F("MC_DecayGammaEtap_Pt","MC_DecayGammaEtap_Pt",250,0,25);
 855                         fMCList[iCut]->Add(fHistoMCDecayGammaEtapPt[iCut]);
 856                         fHistoMCDecayGammaPhiPt[iCut] = new TH1F("MC_DecayGammaPhi_Pt","MC_DecayGammaPhi_Pt",250,0,25);
 857                         fMCList[iCut]->Add(fHistoMCDecayGammaPhiPt[iCut]);
 858                         fHistoMCDecayGammaSigmaPt[iCut] = new TH1F("MC_DecayGammaSigma_Pt","MC_DecayGammaSigma_Pt",250,0,25);
 859                         fMCList[iCut]->Add(fHistoMCDecayGammaSigmaPt[iCut]);
 860                         fHistoMCConvGammaPt[iCut] = new TH1F("MC_ConvGamma_Pt","MC_ConvGamma_Pt",250,0,25);
 861                         fMCList[iCut]->Add(fHistoMCConvGammaPt[iCut]);
 862
 863                         if (fDoPhotonQA > 0){
 864                                 fHistoMCConvGammaR[iCut] = new TH1F("MC_ConvGamma_R","MC_ConvGamma_R",800,0,200);
 865                                 fMCList[iCut]->Add(fHistoMCConvGammaR[iCut]);
 866                                 fHistoMCConvGammaEta[iCut] = new TH1F("MC_ConvGamma_Eta","MC_ConvGamma_Eta",2000,-2,2);
 867                                 fMCList[iCut]->Add(fHistoMCConvGammaEta[iCut]);
 868                         }
 869
 870                         if(fDoMesonAnalysis){
 871                                 fHistoMCPi0Pt[iCut] = new TH1F("MC_Pi0_Pt","MC_Pi0_Pt",250,0,25);
 872                                 fHistoMCPi0Pt[iCut]->Sumw2();
 873                                 fMCList[iCut]->Add(fHistoMCPi0Pt[iCut]);
 874                                 fHistoMCPi0WOWeightPt[iCut] = new TH1F("MC_Pi0_WOWeights_Pt","MC_Pi0_WOWeights_Pt",250,0,25);
 875                                 fHistoMCPi0WOWeightPt[iCut]->Sumw2();
 876                                 fMCList[iCut]->Add(fHistoMCPi0WOWeightPt[iCut]);
 877
 878                                 fHistoMCEtaPt[iCut] = new TH1F("MC_Eta_Pt","MC_Eta_Pt",250,0,25);
 879                                 fHistoMCEtaPt[iCut]->Sumw2();
 880                                 fMCList[iCut]->Add(fHistoMCEtaPt[iCut]);
 881                                 fHistoMCEtaWOWeightPt[iCut] = new TH1F("MC_Eta_WOWeights_Pt","MC_Eta_WOWeights_Pt",250,0,25);
 882                                 fHistoMCEtaWOWeightPt[iCut]->Sumw2();
 883                                 fMCList[iCut]->Add(fHistoMCEtaWOWeightPt[iCut]);
 884                                 fHistoMCPi0InAccPt[iCut] = new TH1F("MC_Pi0InAcc_Pt","MC_Pi0InAcc_Pt",250,0,25);
 885                                 fHistoMCPi0InAccPt[iCut]->Sumw2();
 886                                 fMCList[iCut]->Add(fHistoMCPi0InAccPt[iCut]);
 887                                 fHistoMCEtaInAccPt[iCut] = new TH1F("MC_EtaInAcc_Pt","MC_EtaInAcc_Pt",250,0,25);
 888                                 fHistoMCEtaInAccPt[iCut]->Sumw2();
 889                                 fMCList[iCut]->Add(fHistoMCEtaInAccPt[iCut]);
 890                                 if (fDoMesonQA > 0){
 891                                         fHistoMCPi0PtY[iCut] = new TH2F("MC_Pi0_Pt_Y","MC_Pi0_Pt_Y",150,0.03,15.,150,-1.5,1.5);
 892                                         fHistoMCPi0PtY[iCut]->Sumw2();
 893                                         SetLogBinningXTH2(fHistoMCPi0PtY[iCut]);
 894                                         fMCList[iCut]->Add(fHistoMCPi0PtY[iCut]);
 895                                         fHistoMCEtaPtY[iCut] = new TH2F("MC_Eta_Pt_Y","MC_Eta_Pt_Y",150,0.03,15.,150,-1.5,1.5);
 896                                         fHistoMCEtaPtY[iCut]->Sumw2();
 897                                         SetLogBinningXTH2(fHistoMCEtaPtY[iCut]);
 898                                         fMCList[iCut]->Add(fHistoMCEtaPtY[iCut]);
 899                                         fHistoMCK0sPt[iCut] = new TH1F("MC_K0s_Pt","MC_K0s_Pt",150,0,15);
 900                                         fHistoMCK0sPt[iCut]->Sumw2();
 901                                         fMCList[iCut]->Add(fHistoMCK0sPt[iCut]);
 902                                         fHistoMCK0sWOWeightPt[iCut] = new TH1F("MC_K0s_WOWeights_Pt","MC_K0s_WOWeights_Pt",150,0,15);
 903                                         fHistoMCK0sWOWeightPt[iCut]->Sumw2();
 904                                         fMCList[iCut]->Add(fHistoMCK0sWOWeightPt[iCut]);
 905                                         fHistoMCK0sPtY[iCut] = new TH2F("MC_K0s_Pt_Y","MC_K0s_Pt_Y",150,0.03,15.,150,-1.5,1.5);
 906                                         fHistoMCK0sPtY[iCut]->Sumw2();
 907                                         SetLogBinningXTH2(fHistoMCK0sPtY[iCut]);
 908                                         fMCList[iCut]->Add(fHistoMCK0sPtY[iCut]);
 909
 910                                         fHistoMCSecPi0Source[iCut] = new TH1F("MC_SecPi0_Source","MC_SecPi0_Source",5000,0.,5000);
 911                                         fMCList[iCut]->Add(fHistoMCSecPi0Source[iCut]);
 912                                         fHistoMCSecEtaSource[iCut] = new TH1F("MC_SecEta_Source","MC_SecEta_Source",5000,0,5000);
 913                                         fMCList[iCut]->Add(fHistoMCSecEtaSource[iCut]);
 914                                         fHistoMCSecPi0PtvsSource[iCut] = new TH2F("MC_SecPi0_Pt_Source","MC_SecPi0_Pt_Source",250,0.0,25.,16,-0.5,15.5);
 915                                         fHistoMCSecPi0PtvsSource[iCut]->Sumw2();
 916                                         fMCList[iCut]->Add(fHistoMCSecPi0PtvsSource[iCut]);
 917                                         fHistoMCSecEtaPt[iCut] = new TH1F("MC_SecEta_Pt","MC_SecEta_Pt",250,0,25);
 918                                         fHistoMCSecEtaPt[iCut]->Sumw2();
 919                                         fMCList[iCut]->Add(fHistoMCSecEtaPt[iCut]);
 920                                 }
 921
 922                         }
 923                         fTrueList[iCut] = new TList();
 924                         fTrueList[iCut]->SetName(Form("%s_%s_%s_%s True histograms",cutstringEvent.Data(),cutstringPhoton.Data(),cutstringCalo.Data(),cutstringMeson.Data()));
 925                         fTrueList[iCut]->SetOwner(kTRUE);
 926                         fCutFolder[iCut]->Add(fTrueList[iCut]);
 927
 928                         fHistoTrueConvGammaPt[iCut] = new TH1F("ESD_TrueConvGamma_Pt","ESD_TrueConvGamma_Pt",250,0,25);
 929                         fTrueList[iCut]->Add(fHistoTrueConvGammaPt[iCut]);
 930
 931                         fHistoTrueConvPi0GammaPt[iCut] = new TH1F("ESD_TrueConvGamma_Pt","ESD_TrueConvGamma_Pt",250,0,25);
 932                         fTrueList[iCut]->Add(fHistoTrueConvPi0GammaPt[iCut]);
 933
 934                         fHistoCombinatorialPt[iCut] = new TH2F("ESD_TrueCombinatorial_Pt","ESD_TrueCombinatorial_Pt",250,0,25,16,-0.5,15.5);
 935                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 1,"Elec+Elec");
 936                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 2,"Elec+Pion");
 937                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 3,"Elec+Kaon");
 938                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 4,"Elec+Proton");
 939                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 5,"Elec+Muon");
 940                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 6,"Pion+Pion");
 941                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 7,"Pion+Kaon");
 942                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 8,"Pion+Proton");
 943                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 9,"Pion+Muon");
 944                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(10,"Kaon+Kaon");
 945                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(11,"Kaon+Proton");
 946                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(12,"Kaon+Muon");
 947                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(13,"Proton+Proton");
 948                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(14,"Proton+Muon");
 949                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(15,"Muon+Muon");
 950                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(16,"Rest");
 951                         fTrueList[iCut]->Add(fHistoCombinatorialPt[iCut]);
 952                         fHistoTruePrimaryConvGammaPt[iCut] = new TH1F("ESD_TruePrimaryConvGamma_Pt","ESD_TruePrimaryConvGamma_Pt",250,0,25);
 953                         fTrueList[iCut]->Add(fHistoTruePrimaryConvGammaPt[iCut]);
 954                         fHistoTrueSecondaryConvGammaPt[iCut] = new TH1F("ESD_TrueSecondaryConvGamma_Pt","ESD_TrueSecondaryConvGamma_Pt",250,0,25);
 955                         fTrueList[iCut]->Add(fHistoTrueSecondaryConvGammaPt[iCut]);
 956
 957                         fHistoTrueSecondaryConvGammaFromXFromK0sPt[iCut] = new TH1F("ESD_TrueSecondaryConvGammaFromXFromK0s_Pt", "ESD_TrueSecondaryConvGammaFromXFromK0s_Pt",250,0,25);
 958                         fTrueList[iCut]->Add(fHistoTrueSecondaryConvGammaFromXFromK0sPt[iCut]);
 959                         fHistoTrueSecondaryConvGammaFromXFromLambdaPt[iCut] = new TH1F("ESD_TrueSecondaryConvGammaFromXFromLambda_Pt", "ESD_TrueSecondaryConvGammaFromXFromLambda_Pt",250,0,25);
 960                         fTrueList[iCut]->Add(fHistoTrueSecondaryConvGammaFromXFromLambdaPt[iCut]);
 961
 962                         fHistoTrueDalitzPsiPairDeltaPhi[iCut] = new TH2F("ESD_TrueDalitzPsiPairDeltaPhi_Pt", "ESD_TrueDalitzPsiPairDeltaPhi_Pt",100,-0.5,2,100,-0.5,0.5);
 963                         fTrueList[iCut]->Add(fHistoTrueDalitzPsiPairDeltaPhi[iCut]);
 964
 965                         fHistoTrueGammaPsiPairDeltaPhi[iCut] = new TH2F("ESD_TrueGammaPsiPairDeltaPhi_Pt", "ESD_TrueGammaPsiPairDeltaPhi_Pt",100,-0.5,2,100,-0.5,0.5);
 966                         fTrueList[iCut]->Add(fHistoTrueGammaPsiPairDeltaPhi[iCut]);
 967
 968                         fHistoTruePrimaryConvGammaESDPtMCPt[iCut] = new TH2F("ESD_TruePrimaryConvGammaESD_PtMCPt", "ESD_TruePrimaryConvGammaESD_PtMCPt",250,0,25,250,0,25);
 969                         fTrueList[iCut]->Add(fHistoTruePrimaryConvGammaESDPtMCPt[iCut]);
 970
 971                         fHistoTrueClusGammaPt[iCut] = new TH1F("TrueClusGamma_Pt","ESD_TruePrimaryClusGamma_Pt",250,0,25);
 972                         fTagOutputList[iCut]->Add(fHistoTrueClusGammaPt[iCut]);
 973                         fHistoTruePrimaryClusGammaPt[iCut] = new TH1F("TruePrimaryClusGamma_Pt","ESD_TruePrimaryClusGamma_Pt",250,0,25);
 974                         fTagOutputList[iCut]->Add(fHistoTruePrimaryClusGammaPt[iCut]);
 975                         fHistoTruePrimaryClusGammaESDPtMCPt[iCut] = new TH2F("TruePrimaryClusGamma_Pt_MCPt","ESD_TruePrimaryClusGamma_MCPt",250,0,25,250,0,25);
 976                         fTagOutputList[iCut]->Add(fHistoTruePrimaryClusGammaESDPtMCPt[iCut]);
 977
 978                         if(fDoMesonAnalysis){
 979                                 fHistoTrueMotherInvMassPt[iCut] = new TH2F("ESD_TrueMother_InvMass_Pt","ESD_TrueMother_InvMass_Pt",800,0,0.8,250,0,25);
 980                                 fTrueList[iCut]->Add(fHistoTrueMotherInvMassPt[iCut]);
 981                                 fHistoTruePrimaryMotherInvMassPt[iCut] = new TH2F("ESD_TruePrimaryMother_InvMass_Pt", "ESD_TruePrimaryMother_InvMass_Pt", 800,0,0.8,250,0,25);
 982                                 fHistoTruePrimaryMotherInvMassPt[iCut]->Sumw2();
 983                                 fTrueList[iCut]->Add(fHistoTruePrimaryMotherInvMassPt[iCut]);
 984                                 fHistoTruePrimaryMotherW0WeightingInvMassPt[iCut] = new TH2F("ESD_TruePrimaryMotherW0Weights_InvMass_Pt", "ESD_TruePrimaryMotherW0Weights_InvMass_Pt", 800,0,0.8,250,0,25);
 985                                 fHistoTruePrimaryMotherW0WeightingInvMassPt[iCut]->Sumw2();
 986                                 fTrueList[iCut]->Add(fHistoTruePrimaryMotherW0WeightingInvMassPt[iCut]);
 987                                 fProfileTruePrimaryMotherWeightsInvMassPt[iCut] = new TProfile2D("ESD_TruePrimaryMotherWeights_InvMass_Pt", "ESD_TruePrimaryMotherWeights_InvMass_Pt", 800,0,0.8,250,0,25);
 988                                 fProfileTruePrimaryMotherWeightsInvMassPt[iCut]->Sumw2();
 989                                 fTrueList[iCut]->Add(fProfileTruePrimaryMotherWeightsInvMassPt[iCut]);
 990                                 fHistoTrueSecondaryMotherInvMassPt[iCut] = new TH2F("ESD_TrueSecondaryMother_InvMass_Pt", "ESD_TrueSecondaryMother_InvMass_Pt", 800,0,0.8,250,0,25);
 991                                 fHistoTrueSecondaryMotherInvMassPt[iCut]->Sumw2();
 992                                 fTrueList[iCut]->Add(fHistoTrueSecondaryMotherInvMassPt[iCut]);
 993                                 fHistoTrueSecondaryMotherFromK0sInvMassPt[iCut] = new TH2F("ESD_TrueSecondaryMotherFromK0s_InvMass_Pt","ESD_TrueSecondaryMotherFromK0s_InvMass_Pt",800,0,0.8,250,0,25);
 994                                 fHistoTrueSecondaryMotherFromK0sInvMassPt[iCut]->Sumw2();
 995                                 fTrueList[iCut]->Add(fHistoTrueSecondaryMotherFromK0sInvMassPt[iCut]);
 996                                 fHistoTrueSecondaryMotherFromEtaInvMassPt[iCut] = new TH2F("ESD_TrueSecondaryMotherFromEta_InvMass_Pt","ESD_TrueSecondaryMotherFromEta_InvMass_Pt",800,0,0.8,250,0,25);
 997                                 fTrueList[iCut]->Add(fHistoTrueSecondaryMotherFromEtaInvMassPt[iCut]);
 998                                 fHistoTrueSecondaryMotherFromLambdaInvMassPt[iCut] = new TH2F("ESD_TrueSecondaryMotherFromLambda_InvMass_Pt","ESD_TrueSecondaryMotherFromLambda_InvMass_Pt",800,0,0.8,250,0,25);
 999                                 fTrueList[iCut]->Add(fHistoTrueSecondaryMotherFromLambdaInvMassPt[iCut]);
1000                                 fHistoTrueMotherDalitzInvMassPt[iCut] = new TH2F("ESD_TrueDalitz_InvMass_Pt","ESD_TrueDalitz_InvMass_Pt",800,0,0.8,250,0,25);
1001                                 fTrueList[iCut]->Add(fHistoTrueMotherDalitzInvMassPt[iCut]);
1002                                 if (fDoMesonQA > 0){
1003                                         fHistoTruePrimaryPi0MCPtResolPt[iCut] = new TH2F("ESD_TruePrimaryPi0_MCPt_ResolPt","ESD_TruePrimaryPi0_ResolPt_MCPt",500,0.03,25,1000,-1.,1.);
1004                                         fHistoTruePrimaryPi0MCPtResolPt[iCut]->Sumw2();
1005                                         SetLogBinningXTH2(fHistoTruePrimaryPi0MCPtResolPt[iCut]);
1006                                         fTrueList[iCut]->Add(fHistoTruePrimaryPi0MCPtResolPt[iCut]);
1007                                         fHistoTruePrimaryEtaMCPtResolPt[iCut]  = new TH2F("ESD_TruePrimaryEta_MCPt_ResolPt","ESD_TruePrimaryEta_ResolPt_MCPt",500,0.03,25,1000,-1.,1.);
1008                                         fHistoTruePrimaryEtaMCPtResolPt[iCut]->Sumw2();
1009                                         SetLogBinningXTH2(fHistoTruePrimaryEtaMCPtResolPt[iCut]);
1010                                         fTrueList[iCut]->Add(fHistoTruePrimaryEtaMCPtResolPt[iCut]);
1011                                         fHistoTrueBckGGInvMassPt[iCut] = new TH2F("ESD_TrueBckGG_InvMass_Pt","ESD_TrueBckGG_InvMass_Pt",800,0,0.8,250,0,25);
1012                                         fTrueList[iCut]->Add(fHistoTrueBckGGInvMassPt[iCut]);
1013                                         fHistoTrueBckContInvMassPt[iCut] = new TH2F("ESD_TrueBckCont_InvMass_Pt","ESD_TrueBckCont_InvMass_Pt",800,0,0.8,250,0,25);
1014                                         fTrueList[iCut]->Add(fHistoTrueBckContInvMassPt[iCut]);
1015                                         fHistoTrueK0sWithPi0DaughterMCPt[iCut] = new TH1F("ESD_TrueK0sWithPi0Daughter_MCPt","ESD_TrueK0sWithPi0Daughter_MCPt",250,0,25);
1016                                         fTrueList[iCut]->Add(fHistoTrueK0sWithPi0DaughterMCPt[iCut]);
1017                                         fHistoTrueEtaWithPi0DaughterMCPt[iCut] = new TH1F("ESD_TrueEtaWithPi0Daughter_MCPt","ESD_TrueEtaWithPi0Daughter_MCPt",250,0,25);
1018                                         fTrueList[iCut]->Add(fHistoTrueEtaWithPi0DaughterMCPt[iCut]);
1019                                         fHistoTrueLambdaWithPi0DaughterMCPt[iCut] = new TH1F("ESD_TrueLambdaWithPi0Daughter_MCPt","ESD_TrueLambdaWithPi0Daughter_MCPt",250,0,25);
1020                                         fTrueList[iCut]->Add(fHistoTrueLambdaWithPi0DaughterMCPt[iCut]);
1021
1022                                         fHistoTruePi0PtY[iCut] = new TH2F("ESD_TruePi0_Pt_Y","ESD_TruePi0_Pt_Y",150,0.03,15.,150,-1.5,1.5);
1023                                         SetLogBinningXTH2(fHistoTruePi0PtY[iCut]);
1024                                         fTrueList[iCut]->Add(fHistoTruePi0PtY[iCut]);
1025                                         fHistoTrueEtaPtY[iCut] = new TH2F("ESD_TrueEta_Pt_Y","ESD_TrueEta_Pt_Y",150,0.03,15.,150,-1.5,1.5);
1026                                         SetLogBinningXTH2(fHistoTrueEtaPtY[iCut]);
1027                                         fTrueList[iCut]->Add(fHistoTrueEtaPtY[iCut]);
1028                                         fHistoTruePi0PtAlpha[iCut] = new TH2F("ESD_TruePi0_Pt_Alpha","ESD_TruePi0_Pt_Alpha",150,0.03,15.,100,0,1);
1029                                         SetLogBinningXTH2(fHistoTruePi0PtAlpha[iCut]);
1030                                         fTrueList[iCut]->Add(fHistoTruePi0PtAlpha[iCut]);
1031                                         fHistoTrueEtaPtAlpha[iCut] = new TH2F("ESD_TrueEta_Pt_Alpha","ESD_TrueEta_Pt_Alpha",150,0.03,15.,100,0,1);
1032                                         SetLogBinningXTH2(fHistoTrueEtaPtAlpha[iCut]);
1033                                         fTrueList[iCut]->Add(fHistoTrueEtaPtAlpha[iCut]);
1034
1035                                         fHistoTruePi0PtOpenAngle[iCut] = new TH2F("ESD_TruePi0_Pt_OpenAngle","ESD_TruePi0_Pt_OpenAngle",150,0.03,15.,200,0,2*TMath::Pi());
1036                                         SetLogBinningXTH2(fHistoTruePi0PtOpenAngle[iCut]);
1037                                         fTrueList[iCut]->Add(fHistoTruePi0PtOpenAngle[iCut]);
1038                                         fHistoTrueEtaPtOpenAngle[iCut] = new TH2F("ESD_TrueEta_Pt_OpenAngle","ESD_TrueEta_Pt_OpenAngle",150,0.03,15.,200,0,2*TMath::Pi());
1039                                         SetLogBinningXTH2(fHistoTrueEtaPtOpenAngle[iCut]);
1040                                         fTrueList[iCut]->Add(fHistoTrueEtaPtOpenAngle[iCut]);
1041
1042                                         fHistoTrueConvGammaEta[iCut] = new TH1F("ESD_TrueConvGamma_Eta","ESD_TrueConvGamma_Eta",2000,-2,2);
1043                                         fTrueList[iCut]->Add(fHistoTrueConvGammaEta[iCut]);
1044
1045                                 }
1046                         }
1047                 }
1048         }
1049
1050         fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
1051         if(!fV0Reader){printf("Error: No V0 Reader");return;} // GetV0Reader
1052
1053         if(fV0Reader)
1054                 if((AliConversionPhotonCuts*)fV0Reader->GetConversionCuts())
1055                         if(((AliConversionPhotonCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms())
1056                                 fOutputContainer->Add(((AliConversionPhotonCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms());
1057         if(fV0Reader)
1058                 if((AliConvEventCuts*)fV0Reader->GetEventCuts())
1059                         if(((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetCutHistograms())
1060                                 fOutputContainer->Add(((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetCutHistograms());
1061
1062
1063         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
1064                 if(!((AliConvEventCuts*)fEventCutArray->At(iCut))) continue;
1065                 if(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutHistograms()){
1066                         fCutFolder[iCut]->Add(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutHistograms());
1067                 }
1068                 if(!((AliConversionPhotonCuts*)fCutArray->At(iCut))) continue;
1069                 if(((AliConversionPhotonCuts*)fCutArray->At(iCut))->GetCutHistograms()){
1070                         fCutFolder[iCut]->Add(((AliConversionPhotonCuts*)fCutArray->At(iCut))->GetCutHistograms());
1071                 }
1072                 if(!((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))) continue;
1073                 if(((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutHistograms()){
1074                         fCutFolder[iCut]->Add(((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutHistograms());
1075                 }
1076                 if(fDoMesonAnalysis){
1077                         if(!((AliConversionMesonCuts*)fMesonCutArray->At(iCut))) continue;
1078                         if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutHistograms()){
1079                                 fCutFolder[iCut]->Add(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutHistograms());
1080                         }
1081                 }
1082         }
1083         PostData(1, fOutputContainer);
1084 }
1085 //_____________________________________________________________________________
1086 Bool_t AliAnalysisTaskGammaConvCalo::Notify()
1087 {
1088         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
1089                 if(!((AliConvEventCuts*)fEventCutArray->At(iCut))->GetDoEtaShift()){
1090                         fProfileEtaShift[iCut]->Fill(0.,(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift()));
1091                         continue; // No Eta Shift requested, continue
1092                 }
1093                 if(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift() == 0.0){ // Eta Shift requested but not set, get shift automatically
1094                         ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCorrectEtaShiftFromPeriod(fV0Reader->GetPeriodName());
1095                         fProfileEtaShift[iCut]->Fill(0.,(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift()));
1096                         ((AliConvEventCuts*)fEventCutArray->At(iCut))->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
1097                         continue;
1098                 }
1099                 else{
1100                         printf(" Gamma Conversion Task %s :: Eta Shift Manually Set to %f \n\n",
1101                                         (((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber()).Data(),((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift());
1102                         fProfileEtaShift[iCut]->Fill(0.,(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift()));
1103                         ((AliConvEventCuts*)fEventCutArray->At(iCut))->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
1104                 }
1105         }
1106
1107         return kTRUE;
1108 }
1109 //_____________________________________________________________________________
1110 void AliAnalysisTaskGammaConvCalo::UserExec(Option_t *)
1111 {
1112         //
1113         // Called for each event
1114         //
1115         Int_t eventQuality = ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEventQuality();
1116         if(eventQuality == 2 || eventQuality == 3){// Event Not Accepted due to MC event missing or wrong trigger for V0ReaderV1
1117                 for(Int_t iCut = 0; iCut<fnCuts; iCut++){
1118                 fHistoNEvents[iCut]->Fill(eventQuality);
1119                 }
1120                 return;
1121         }
1122
1123         if(fIsMC) fMCEvent = MCEvent();
1124         if(fMCEvent == NULL) fIsMC = kFALSE;
1125
1126         fInputEvent = InputEvent();
1127
1128         if(fIsMC && fInputEvent->IsA()==AliESDEvent::Class()){
1129                 fMCStack = fMCEvent->Stack();
1130                 if(fMCStack == NULL) fIsMC = kFALSE;
1131         }
1132
1133         fReaderGammas = fV0Reader->GetReconstructedGammas(); // Gammas from default Cut
1134
1135         // ------------------- BeginEvent ----------------------------
1136
1137         AliEventplane *EventPlane = fInputEvent->GetEventplane();
1138         if(fIsHeavyIon ==1)fEventPlaneAngle = EventPlane->GetEventplane("V0",fInputEvent,2);
1139         else fEventPlaneAngle=0.0;
1140
1141         if(fIsMC && fInputEvent->IsA()==AliAODEvent::Class() && !(fV0Reader->AreAODsRelabeled())){
1142                 RelabelAODPhotonCandidates(kTRUE);    // In case of AODMC relabeling MC
1143                 fV0Reader->RelabelAODs(kTRUE);
1144         }
1145
1146
1147         for(Int_t iCut = 0; iCut<fnCuts; iCut++){
1148
1149                 fiCut = iCut;
1150                 Int_t eventNotAccepted = ((AliConvEventCuts*)fEventCutArray->At(iCut))->IsEventAcceptedByCut(fV0Reader->GetEventCuts(),fInputEvent,fMCEvent,fIsHeavyIon);
1151
1152                 if(eventNotAccepted){
1153                 // cout << "event rejected due to wrong trigger: " <<eventNotAccepted << endl;
1154                         fHistoNEvents[iCut]->Fill(eventNotAccepted); // Check Centrality, PileUp, SDD and V0AND --> Not Accepted => eventQuality = 1
1155                         continue;
1156                 }
1157
1158                 if(eventQuality != 0){// Event Not Accepted
1159                         //cout << "event rejected due to: " <<eventQuality << endl;
1160                         fHistoNEvents[iCut]->Fill(eventQuality);
1161                         continue;
1162                 }
1163
1164                 fHistoNEvents[iCut]->Fill(eventQuality); // Should be 0 here
1165                 fHistoNGoodESDTracks[iCut]->Fill(fV0Reader->GetNumberOfPrimaryTracks());
1166                 if(((AliConvEventCuts*)fEventCutArray->At(iCut))->IsHeavyIon() == 2)    fHistoNV0Tracks[iCut]->Fill(fInputEvent->GetVZEROData()->GetMTotV0A());
1167                         else fHistoNV0Tracks[iCut]->Fill(fInputEvent->GetVZEROData()->GetMTotV0A()+fInputEvent->GetVZEROData()->GetMTotV0C());
1168
1169                 if(fIsMC){
1170                         // Process MC Particle
1171                         if(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetSignalRejection() != 0){
1172                                 if(fInputEvent->IsA()==AliESDEvent::Class()){
1173                                 ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetNotRejectedParticles(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetSignalRejection(),
1174                                                                                                                                                                         ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetAcceptedHeader(),
1175                                                                                                                                                                         fMCEvent);
1176                                 }
1177                                 else if(fInputEvent->IsA()==AliAODEvent::Class()){
1178                                 ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetNotRejectedParticles(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetSignalRejection(),
1179                                                                                                                                                                         ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetAcceptedHeader(),
1180                                                                                                                                                                         fInputEvent);
1181                                 }
1182
1183                                 if(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetAcceptedHeader()){
1184                                         for(Int_t i = 0;i<(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetAcceptedHeader())->GetEntries();i++){
1185                                                 TString nameBin= fHistoMCHeaders[iCut]->GetXaxis()->GetBinLabel(i+1);
1186                                                 if (nameBin.CompareTo("")== 0){
1187                                                         TString nameHeader = ((TObjString*)((TList*)((AliConvEventCuts*)fEventCutArray->At(iCut))
1188                                                                                                                                 ->GetAcceptedHeader())->At(i))->GetString();
1189                                                         fHistoMCHeaders[iCut]->GetXaxis()->SetBinLabel(i+1,nameHeader.Data());
1190                                                 }
1191                                         }
1192                                 }
1193                         }
1194                 }
1195                 if(fIsMC){
1196                 if(fInputEvent->IsA()==AliESDEvent::Class())
1197                         ProcessMCParticles();
1198                 if(fInputEvent->IsA()==AliAODEvent::Class())
1199                         ProcessAODMCParticles();
1200                 }
1201
1202                 // 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)
1203                 ProcessClusters();                                      // process calo clusters
1204                 ProcessPhotonCandidates();                              // Process this cuts gammas
1205
1206                 fHistoNGammaCandidates[iCut]->Fill(fGammaCandidates->GetEntries());
1207                 fHistoNGoodESDTracksVsNGammaCanditates[iCut]->Fill(fV0Reader->GetNumberOfPrimaryTracks(),fGammaCandidates->GetEntries());
1208                 if(fDoMesonAnalysis){ // Meson Analysis
1209                 if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseMCPSmearing() && fIsMC){
1210                         fUnsmearedPx = new Double_t[fGammaCandidates->GetEntries()]; // Store unsmeared Momenta
1211                         fUnsmearedPy = new Double_t[fGammaCandidates->GetEntries()];
1212                         fUnsmearedPz = new Double_t[fGammaCandidates->GetEntries()];
1213                         fUnsmearedE =  new Double_t[fGammaCandidates->GetEntries()];
1214
1215                         for(Int_t gamma=0;gamma<fGammaCandidates->GetEntries();gamma++){ // Smear the AODPhotons in MC
1216                         fUnsmearedPx[gamma] = ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->Px();
1217                         fUnsmearedPy[gamma] = ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->Py();
1218                         fUnsmearedPz[gamma] = ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->Pz();
1219                         fUnsmearedE[gamma] =  ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->E();
1220                         ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->SmearParticle(dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(gamma)));
1221                         }
1222                 }
1223
1224                 PhotonTagging(); // tag PCM photons with calorimeter
1225
1226                 CalculatePi0Candidates(); // Combine Gammas from conversion and from calo
1227
1228                 if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->DoBGCalculation()){
1229                         if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->BackgroundHandlerType() == 0){
1230                         CalculateBackground(); // Combinatorial Background
1231                         UpdateEventByEventData(); // Store Event for mixed Events
1232                         }
1233                         else{
1234                         CalculateBackgroundRP(); // Combinatorial Background
1235                         fBGHandlerRP[iCut]->AddEvent(fGammaCandidates,fInputEvent); // Store Event for mixed Events
1236                         fBGClusHandlerRP[iCut]->AddEvent(fClusterCandidates,fInputEvent); // Store Event for mixed Events
1237                         }
1238                 }
1239                 if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseMCPSmearing() && fIsMC){
1240                         for(Int_t gamma=0;gamma<fGammaCandidates->GetEntries();gamma++){ // Smear the AODPhotons in MC
1241                         ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetPx(fUnsmearedPx[gamma]); // Reset Unsmeared Momenta
1242                         ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetPy(fUnsmearedPy[gamma]);
1243                         ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetPz(fUnsmearedPz[gamma]);
1244                         ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetE(fUnsmearedE[gamma]);
1245                         }
1246                         delete[] fUnsmearedPx; fUnsmearedPx = 0x0;
1247                         delete[] fUnsmearedPy; fUnsmearedPy = 0x0;
1248                         delete[] fUnsmearedPz; fUnsmearedPz = 0x0;
1249                         delete[] fUnsmearedE;  fUnsmearedE  = 0x0;
1250                 }
1251                 }
1252
1253                 fGammaCandidates->Clear(); // delete this cuts good gammas
1254                 fClusterCandidates->Clear(); // delete cluster candidates
1255         }
1256
1257         if(fIsMC && fInputEvent->IsA()==AliAODEvent::Class() && !(fV0Reader->AreAODsRelabeled())){
1258                 RelabelAODPhotonCandidates(kFALSE); // Back to ESDMC Label
1259                 fV0Reader->RelabelAODs(kFALSE);
1260         }
1261
1262         PostData(1, fOutputContainer);
1263 }
1264
1265 //________________________________________________________________________
1266 void AliAnalysisTaskGammaConvCalo::ProcessClusters()
1267 {
1268
1269         Int_t nclus = 0;
1270         nclus = fInputEvent->GetNumberOfCaloClusters();
1271
1272 //      cout << nclus << endl;
1273
1274         if(nclus == 0)  return;
1275
1276         // vertex
1277         Double_t vertex[3] = {0};
1278         InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1279
1280         // Loop over EMCal clusters
1281         for(Int_t i = 0; i < nclus; i++){
1282
1283                 AliVCluster* clus = fInputEvent->GetCaloCluster(i);
1284                 if (!clus) continue;
1285                 if(!((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->ClusterIsSelected(clus,fInputEvent,fIsMC)) continue;
1286                 // TLorentzvector with cluster
1287                 TLorentzVector clusterVector;
1288                 clus->GetMomentum(clusterVector,vertex);
1289
1290                 TLorentzVector* tmpvec = new TLorentzVector();
1291                 tmpvec->SetPxPyPzE(clusterVector.Px(),clusterVector.Py(),clusterVector.Pz(),clusterVector.E());
1292
1293                 // convert to AODConversionPhoton
1294                 AliAODConversionPhoton *PhotonCandidate=new AliAODConversionPhoton(tmpvec);
1295                 if(!PhotonCandidate) continue;
1296
1297                 // get MC label
1298                 Int_t mclab[2] = {0,0};
1299                 if(fIsMC){
1300                         mclab[0] = clus->GetLabel();
1301                         mclab[1] = clus->GetLabel();
1302 //                      cout << "mclabels: " << mclab[0] << " " << mclab[1] << endl;
1303                 }
1304
1305                 PhotonCandidate->SetMCLabel(mclab);
1306
1307         //    fIsFromMBHeader = kTRUE;
1308         //    if(fIsMC){
1309         //      Int_t isPosFromMBHeader
1310         //      = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(mclab, fMCStack, fInputEvent);
1311         //      if(isPosFromMBHeader == 0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1312         //
1313         //      Int_t isNegFromMBHeader
1314         //      = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(mclab, fMCStack, fInputEvent);
1315         //      if(isNegFromMBHeader == 0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1316         //
1317         //      if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE;
1318         //    }
1319
1320                 fHistoClusGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
1321                 fClusterCandidates->Add(PhotonCandidate); // if no second loop is required add to events good gammas
1322
1323         //    if(fIsFromMBHeader){
1324         //      fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
1325         //      if (fDoPhotonQA > 0){
1326         //        fHistoConvGammaR[fiCut]->Fill(PhotonCandidate->GetConversionRadius());
1327         //        fHistoConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
1328         //      }
1329         //    }
1330
1331                 if(fIsMC){
1332                         ProcessTrueClusterCandidates(PhotonCandidate);
1333                 }
1334
1335         //    if (fIsFromMBHeader && fDoPhotonQA == 2){
1336         //      if (fIsHeavyIon == 1 && PhotonCandidate->Pt() > 0.399 && PhotonCandidate->Pt() < 12.){
1337         //        fPtGamma = PhotonCandidate->Pt();
1338         //        fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1339         //        fRConvPhoton = PhotonCandidate->GetConversionRadius();
1340         //        fEtaPhoton = PhotonCandidate->GetPhotonEta();
1341         //        fCharCatPhoton = PhotonCandidate->GetPhotonQuality();
1342         //        fTreeConvGammaPtDcazCat[fiCut]->Fill();
1343         //    } else if ( PhotonCandidate->Pt() > 0.299 && PhotonCandidate->Pt() < 16.){
1344         //        fPtGamma = PhotonCandidate->Pt();
1345         //        fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1346         //        fRConvPhoton = PhotonCandidate->GetConversionRadius();
1347         //        fEtaPhoton = PhotonCandidate->GetPhotonEta();
1348         //        fCharCatPhoton = PhotonCandidate->GetPhotonQuality();
1349         //        fTreeConvGammaPtDcazCat[fiCut]->Fill();
1350         //      }
1351         //    }
1352                 delete tmpvec;
1353         }
1354
1355 }
1356
1357 //________________________________________________________________________
1358 void AliAnalysisTaskGammaConvCalo::ProcessTrueClusterCandidates(AliAODConversionPhoton *TruePhotonCandidate)
1359 {
1360
1361         TParticle *Photon = TruePhotonCandidate->GetMCParticle(fMCStack);
1362
1363         if(Photon == NULL){
1364         //    cout << "no photon" << endl;
1365                 return;
1366         }
1367
1368         if(Photon->GetPdgCode() != 22){
1369         //    cout << "pdg code: " << Photon->GetPdgCode() << endl;
1370                 return; // Particle is no Photon
1371         }
1372
1373         // True Photon
1374         if(fIsFromMBHeader){
1375                 fHistoTrueClusGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1376                 if (fDoPhotonQA > 0) fHistoTrueConvGammaEta[fiCut]->Fill(TruePhotonCandidate->Eta());
1377         }
1378
1379         if(Photon->GetMother(0) <= fMCStack->GetNprimary()){
1380                 // Count just primary MC Gammas as true --> For Ratio esdtruegamma / mcconvgamma
1381                 if(fIsFromMBHeader){
1382                         fHistoTruePrimaryClusGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1383                         fHistoTruePrimaryClusGammaESDPtMCPt[fiCut]->Fill(TruePhotonCandidate->Pt(),Photon->Pt()); // Allways Filled
1384                 }
1385         }
1386
1387         // maybe later
1388         //  else{
1389         //    if(fIsFromMBHeader){
1390         //      fCharPhotonMCInfo = 2;
1391         //      fHistoTrueSecondaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1392         //      if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 &&
1393         //         fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 3122){
1394         //        fHistoTrueSecondaryConvGammaFromXFromLambdaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1395         //        fCharPhotonMCInfo = 5;
1396         //      }
1397         //      if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 &&
1398         //         fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 310){
1399         //        fHistoTrueSecondaryConvGammaFromXFromK0sPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1400         //        fCharPhotonMCInfo = 4;
1401         //      }
1402         //      if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 &&
1403         //         fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 221){
1404         //        fCharPhotonMCInfo = 3;
1405         //      }
1406         //    }
1407         //  }
1408
1409 }
1410
1411
1412 //________________________________________________________________________
1413 void AliAnalysisTaskGammaConvCalo::ProcessPhotonCandidates()
1414 {
1415         Int_t nV0 = 0;
1416         TList *GammaCandidatesStepOne = new TList();
1417         TList *GammaCandidatesStepTwo = new TList();
1418         // Loop over Photon Candidates allocated by ReaderV1
1419         for(Int_t i = 0; i < fReaderGammas->GetEntriesFast(); i++){
1420                 AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) fReaderGammas->At(i);
1421                 if(!PhotonCandidate) continue;
1422                 fIsFromMBHeader = kTRUE;
1423                 if(fIsMC && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0){
1424                         Int_t isPosFromMBHeader = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack, fInputEvent);
1425                         if(isPosFromMBHeader == 0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1426                         Int_t isNegFromMBHeader = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack, fInputEvent);
1427                         if(isNegFromMBHeader == 0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1428                         if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE;
1429                 }
1430
1431                 if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelected(PhotonCandidate,fInputEvent)) continue;
1432                 if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(PhotonCandidate->GetPhotonPhi(),fEventPlaneAngle)) continue;
1433                 if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseElecSharingCut() &&
1434                 !((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){
1435                         fGammaCandidates->Add(PhotonCandidate); // if no second loop is required add to events good gammas
1436
1437                         if(fIsFromMBHeader){
1438                                 fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
1439                                 if (fDoPhotonQA > 0){
1440                                 fHistoConvGammaR[fiCut]->Fill(PhotonCandidate->GetConversionRadius());
1441                                 fHistoConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
1442                                 }
1443                         }
1444                         if(fIsMC){
1445                                 if(fInputEvent->IsA()==AliESDEvent::Class())
1446                                 ProcessTruePhotonCandidates(PhotonCandidate);
1447                                 if(fInputEvent->IsA()==AliAODEvent::Class())
1448                                 ProcessTruePhotonCandidatesAOD(PhotonCandidate);
1449                         }
1450                         if (fIsFromMBHeader && fDoPhotonQA == 2){
1451                                 if (fIsHeavyIon == 1 && PhotonCandidate->Pt() > 0.399 && PhotonCandidate->Pt() < 12.){
1452                                         fPtGamma = PhotonCandidate->Pt();
1453                                         fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1454                                         fRConvPhoton = PhotonCandidate->GetConversionRadius();
1455                                         fEtaPhoton = PhotonCandidate->GetPhotonEta();
1456                                         fCharCatPhoton = PhotonCandidate->GetPhotonQuality();
1457                                         fTreeConvGammaPtDcazCat[fiCut]->Fill();
1458                                 } else if ( PhotonCandidate->Pt() > 0.299 && PhotonCandidate->Pt() < 16.){
1459                                         fPtGamma = PhotonCandidate->Pt();
1460                                         fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1461                                         fRConvPhoton = PhotonCandidate->GetConversionRadius();
1462                                         fEtaPhoton = PhotonCandidate->GetPhotonEta();
1463                                         fCharCatPhoton = PhotonCandidate->GetPhotonQuality();
1464                                         fTreeConvGammaPtDcazCat[fiCut]->Fill();
1465                                 }
1466                         }
1467                 } else if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseElecSharingCut()){ // if Shared Electron cut is enabled, Fill array, add to step one
1468                         ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->FillElectonLabelArray(PhotonCandidate,nV0);
1469                         nV0++;
1470                         GammaCandidatesStepOne->Add(PhotonCandidate);
1471                 } else if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseElecSharingCut() &&
1472                                 ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){ // shared electron is disabled, step one not needed -> step two
1473                         GammaCandidatesStepTwo->Add(PhotonCandidate);
1474                 }
1475         }
1476         if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseElecSharingCut()){
1477                 for(Int_t i = 0;i<GammaCandidatesStepOne->GetEntries();i++){
1478                         AliAODConversionPhoton *PhotonCandidate= (AliAODConversionPhoton*) GammaCandidatesStepOne->At(i);
1479                         if(!PhotonCandidate) continue;
1480                         fIsFromMBHeader = kTRUE;
1481                         if(fMCStack && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0){
1482                                 Int_t isPosFromMBHeader = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack, fInputEvent);
1483                                 Int_t isNegFromMBHeader = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack, fInputEvent);
1484                                 if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE;
1485                         }
1486                         if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->RejectSharedElectronV0s(PhotonCandidate,i,GammaCandidatesStepOne->GetEntries())) continue;
1487                         if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){ // To Colse v0s cut diabled, step two not needed
1488                                 fGammaCandidates->Add(PhotonCandidate);
1489                                 if(fIsFromMBHeader){
1490                                         fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
1491                                         if (fDoPhotonQA > 0){
1492                                                 fHistoConvGammaR[fiCut]->Fill(PhotonCandidate->GetConversionRadius());
1493                                                 fHistoConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
1494                                         }
1495                                 }
1496                         }
1497                         if(fIsMC){
1498                                 if(fInputEvent->IsA()==AliESDEvent::Class())
1499                                         ProcessTruePhotonCandidates(PhotonCandidate);
1500                                 if(fInputEvent->IsA()==AliAODEvent::Class())
1501                                         ProcessTruePhotonCandidatesAOD(PhotonCandidate);
1502                         } else GammaCandidatesStepTwo->Add(PhotonCandidate); // Close v0s cut enabled -> add to list two
1503
1504                         if (fIsFromMBHeader && fDoPhotonQA == 2){
1505                                 if (fIsHeavyIon ==1 && PhotonCandidate->Pt() > 0.399 && PhotonCandidate->Pt() < 12.){
1506                                         fPtGamma = PhotonCandidate->Pt();
1507                                         fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1508                                         fRConvPhoton = PhotonCandidate->GetConversionRadius();
1509                                         fEtaPhoton = PhotonCandidate->GetPhotonEta();
1510                                         fCharCatPhoton = PhotonCandidate->GetPhotonQuality();
1511                                         fTreeConvGammaPtDcazCat[fiCut]->Fill();
1512                                 } else if ( PhotonCandidate->Pt() > 0.299 && PhotonCandidate->Pt() < 16.){
1513                                         fPtGamma = PhotonCandidate->Pt();
1514                                         fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1515                                         fRConvPhoton = PhotonCandidate->GetConversionRadius();
1516                                         fEtaPhoton = PhotonCandidate->GetPhotonEta();
1517                                         fCharCatPhoton = PhotonCandidate->GetPhotonQuality();
1518                                         fTreeConvGammaPtDcazCat[fiCut]->Fill();
1519                                 }
1520                         }
1521                 }
1522         }
1523         if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){
1524                 for(Int_t i = 0;i<GammaCandidatesStepTwo->GetEntries();i++){
1525                         AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) GammaCandidatesStepTwo->At(i);
1526                         if(!PhotonCandidate) continue;
1527                         fIsFromMBHeader = kTRUE;
1528                         if(fMCStack && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0){
1529                                 Int_t isPosFromMBHeader = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack, fInputEvent);
1530                                 Int_t isNegFromMBHeader = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack, fInputEvent);
1531                                 if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE;
1532                         }
1533                         if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->RejectToCloseV0s(PhotonCandidate,GammaCandidatesStepTwo,i)) continue;
1534                         fGammaCandidates->Add(PhotonCandidate); // Add gamma to current cut TList
1535                         if(fIsFromMBHeader){
1536                                 fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
1537                                 if (fDoPhotonQA > 0){
1538                                         fHistoConvGammaR[fiCut]->Fill(PhotonCandidate->GetConversionRadius());
1539                                         fHistoConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
1540                                 }
1541                         }
1542                         if(fIsMC){
1543                                 if(fInputEvent->IsA()==AliESDEvent::Class())
1544                                         ProcessTruePhotonCandidates(PhotonCandidate);
1545                                 if(fInputEvent->IsA()==AliAODEvent::Class())
1546                                         ProcessTruePhotonCandidatesAOD(PhotonCandidate);
1547                         }
1548                         if (fIsFromMBHeader){
1549                                 if (fIsHeavyIon == 1 && PhotonCandidate->Pt() > 0.399 && PhotonCandidate->Pt() < 12.){
1550                                         fPtGamma = PhotonCandidate->Pt();
1551                                         fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1552                                         fRConvPhoton = PhotonCandidate->GetConversionRadius();
1553                                         fEtaPhoton = PhotonCandidate->GetPhotonEta();
1554                                         fCharCatPhoton = PhotonCandidate->GetPhotonQuality();
1555                                         fTreeConvGammaPtDcazCat[fiCut]->Fill();
1556                                 } else if ( PhotonCandidate->Pt() > 0.299 && PhotonCandidate->Pt() < 16.){
1557                                         fPtGamma = PhotonCandidate->Pt();
1558                                         fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1559                                         fRConvPhoton = PhotonCandidate->GetConversionRadius();
1560                                         fEtaPhoton = PhotonCandidate->GetPhotonEta();
1561                                         fCharCatPhoton = PhotonCandidate->GetPhotonQuality();
1562                                         fTreeConvGammaPtDcazCat[fiCut]->Fill();
1563                                 }
1564                         }
1565                 }
1566         }
1567
1568         delete GammaCandidatesStepOne;
1569         GammaCandidatesStepOne = 0x0;
1570         delete GammaCandidatesStepTwo;
1571         GammaCandidatesStepTwo = 0x0;
1572
1573 }
1574 //________________________________________________________________________
1575 void AliAnalysisTaskGammaConvCalo::ProcessTruePhotonCandidatesAOD(AliAODConversionPhoton *TruePhotonCandidate)
1576 {
1577
1578         Double_t magField = fInputEvent->GetMagneticField();
1579         if( magField  < 0.0 ){
1580                 magField =  1.0;
1581         }
1582         else {
1583                 magField =  -1.0;
1584         }
1585
1586         TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1587         AliAODMCParticle *posDaughter = (AliAODMCParticle*) AODMCTrackArray->At(TruePhotonCandidate->GetMCLabelPositive());
1588         AliAODMCParticle *negDaughter = (AliAODMCParticle*) AODMCTrackArray->At(TruePhotonCandidate->GetMCLabelNegative());
1589         fCharPhotonMCInfo = 0;
1590
1591         if(posDaughter == NULL || negDaughter == NULL) return; // One particle does not exist
1592         Int_t pdgCode[2] = {abs(posDaughter->GetPdgCode()),abs(negDaughter->GetPdgCode())};
1593
1594         if(posDaughter->GetMother() != negDaughter->GetMother()){
1595                 FillPhotonCombinatorialBackgroundHist(TruePhotonCandidate, pdgCode);
1596                 fCharPhotonMCInfo = 1;
1597                 return;
1598         }
1599         else if(posDaughter->GetMother() == -1){
1600                 FillPhotonCombinatorialBackgroundHist(TruePhotonCandidate, pdgCode);
1601                 fCharPhotonMCInfo = 1;
1602                 return;
1603         }
1604
1605         if(pdgCode[0]!=11 || pdgCode[1]!=11){
1606                 fCharPhotonMCInfo = 1;
1607                 return; //One Particle is not a electron
1608         }
1609         if(posDaughter->GetPdgCode()==negDaughter->GetPdgCode()){
1610                 fCharPhotonMCInfo = 1;
1611                 return; // Same Charge
1612         }
1613
1614         AliAODMCParticle *Photon = (AliAODMCParticle*) AODMCTrackArray->At(posDaughter->GetMother());
1615         AliVTrack * electronCandidate = ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->GetTrack(fInputEvent,TruePhotonCandidate->GetTrackLabelNegative() );
1616         AliVTrack * positronCandidate = ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->GetTrack(fInputEvent,TruePhotonCandidate->GetTrackLabelPositive() );
1617         Double_t deltaPhi = magField * TVector2::Phi_mpi_pi( electronCandidate->Phi()-positronCandidate->Phi());
1618
1619         if(Photon->GetPdgCode() != 22){
1620                 fHistoTrueDalitzPsiPairDeltaPhi[fiCut]->Fill(deltaPhi,TruePhotonCandidate->GetPsiPair());
1621                 fCharPhotonMCInfo = 1;
1622                 return; // Mother is no Photon
1623         }
1624
1625         if(((posDaughter->GetMCProcessCode())) != 5 || ((negDaughter->GetMCProcessCode())) != 5){
1626                 fCharPhotonMCInfo = 1;
1627                 return;// check if the daughters come from a conversion
1628         }
1629         // STILL A BUG IN ALIROOT >>8 HAS TPO BE REMOVED AFTER FIX
1630
1631
1632
1633         // True Photon
1634         if(fIsFromMBHeader){
1635                 fHistoTrueConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1636                 if (fDoPhotonQA > 0) fHistoTrueConvGammaEta[fiCut]->Fill(TruePhotonCandidate->Eta());
1637         }
1638         fHistoTrueGammaPsiPairDeltaPhi[fiCut]->Fill(deltaPhi,TruePhotonCandidate->GetPsiPair());
1639         if(Photon->IsPrimary()){
1640                 // Count just primary MC Gammas as true --> For Ratio esdtruegamma / mcconvgamma
1641                 if(fIsFromMBHeader){
1642                         fCharPhotonMCInfo = 6;
1643                         fHistoTruePrimaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1644                         fHistoTruePrimaryConvGammaESDPtMCPt[fiCut]->Fill(TruePhotonCandidate->Pt(),Photon->Pt()); // Allways Filled
1645                 }
1646                 // (Not Filled for i6, Extra Signal Gamma (parambox) are secondary)
1647         } else {
1648                 if(fIsFromMBHeader){
1649                         fHistoTrueSecondaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1650                         fCharPhotonMCInfo = 2;
1651                         if(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother() > -1 &&
1652                                 ((AliAODMCParticle*)AODMCTrackArray->At(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother()))->GetPdgCode() == 3122){
1653                                 fCharPhotonMCInfo = 5;
1654                                 fHistoTrueSecondaryConvGammaFromXFromLambdaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1655                         }
1656                         if(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother() > -1 &&
1657                                 ((AliAODMCParticle*)AODMCTrackArray->At(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother()))->GetPdgCode() == 310){
1658                                 fCharPhotonMCInfo = 4;
1659                                 fHistoTrueSecondaryConvGammaFromXFromK0sPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1660                         }
1661                         if(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother() > -1 &&
1662                                 ((AliAODMCParticle*)AODMCTrackArray->At(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother()))->GetPdgCode() == 221){
1663                                 fCharPhotonMCInfo = 3;
1664                         }
1665                 }
1666         }
1667
1668 }
1669 //________________________________________________________________________
1670 void AliAnalysisTaskGammaConvCalo::ProcessTruePhotonCandidates(AliAODConversionPhoton *TruePhotonCandidate)
1671 {
1672
1673         Double_t magField = fInputEvent->GetMagneticField();
1674         if( magField  < 0.0 ){
1675                 magField =  1.0;
1676         }
1677         else {
1678                 magField =  -1.0;
1679         }
1680
1681         // Process True Photons
1682         TParticle *posDaughter = TruePhotonCandidate->GetPositiveMCDaughter(fMCStack);
1683         TParticle *negDaughter = TruePhotonCandidate->GetNegativeMCDaughter(fMCStack);
1684         fCharPhotonMCInfo = 0;
1685
1686         if(posDaughter == NULL || negDaughter == NULL) return; // One particle does not exist
1687         Int_t pdgCode[2] = {abs(posDaughter->GetPdgCode()),abs(negDaughter->GetPdgCode())};
1688         fCharPhotonMCInfo = 1;
1689         if(posDaughter->GetMother(0) != negDaughter->GetMother(0)){
1690                 FillPhotonCombinatorialBackgroundHist(TruePhotonCandidate, pdgCode);
1691                 return;
1692         }
1693         else if(posDaughter->GetMother(0) == -1){
1694                 FillPhotonCombinatorialBackgroundHist(TruePhotonCandidate, pdgCode);
1695                 return;
1696         }
1697
1698         if(pdgCode[0]!=11 || pdgCode[1]!=11) return; //One Particle is not a electron
1699
1700         if(posDaughter->GetPdgCode()==negDaughter->GetPdgCode()) return; // Same Charge
1701
1702         TParticle *Photon = TruePhotonCandidate->GetMCParticle(fMCStack);
1703         AliVTrack * electronCandidate = ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->GetTrack(fInputEvent,TruePhotonCandidate->GetTrackLabelNegative() );
1704         AliVTrack * positronCandidate = ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->GetTrack(fInputEvent,TruePhotonCandidate->GetTrackLabelPositive() );
1705         Double_t deltaPhi = magField * TVector2::Phi_mpi_pi( electronCandidate->Phi()-positronCandidate->Phi());
1706
1707         if(Photon->GetPdgCode() != 22){
1708                 fHistoTrueDalitzPsiPairDeltaPhi[fiCut]->Fill(deltaPhi,TruePhotonCandidate->GetPsiPair());
1709                 return; // Mother is no Photon
1710         }
1711
1712         if(posDaughter->GetUniqueID() != 5 || negDaughter->GetUniqueID() !=5) return;// check if the daughters come from a conversion
1713
1714         // True Photon
1715         if(fIsFromMBHeader){
1716                 fHistoTrueConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1717                 if (fDoPhotonQA > 0) fHistoTrueConvGammaEta[fiCut]->Fill(TruePhotonCandidate->Eta());
1718         }
1719         fHistoTrueGammaPsiPairDeltaPhi[fiCut]->Fill(deltaPhi,TruePhotonCandidate->GetPsiPair());
1720         if(posDaughter->GetMother(0) <= fMCStack->GetNprimary()){
1721                 // Count just primary MC Gammas as true --> For Ratio esdtruegamma / mcconvgamma
1722                 if(fIsFromMBHeader){
1723                         fCharPhotonMCInfo = 6;
1724                         fHistoTruePrimaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1725                         fHistoTruePrimaryConvGammaESDPtMCPt[fiCut]->Fill(TruePhotonCandidate->Pt(),Photon->Pt()); // Allways Filled
1726
1727                 }
1728                 // (Not Filled for i6, Extra Signal Gamma (parambox) are secondary)
1729         }
1730         else{
1731                 if(fIsFromMBHeader){
1732                         fCharPhotonMCInfo = 2;
1733                         fHistoTrueSecondaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1734                         if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 &&
1735                                 fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 3122){
1736                                 fHistoTrueSecondaryConvGammaFromXFromLambdaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1737                                 fCharPhotonMCInfo = 5;
1738                         }
1739                         if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 &&
1740                                 fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 310){
1741                                 fHistoTrueSecondaryConvGammaFromXFromK0sPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1742                                 fCharPhotonMCInfo = 4;
1743                         }
1744                         if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 &&
1745                                 fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 221){
1746                                 fCharPhotonMCInfo = 3;
1747                         }
1748                 }
1749         }
1750
1751         // pi0 photon
1752         //Bool_t bpi0 = 0;
1753         Int_t imother = Photon->GetMother(0);
1754         AliMCParticle *McMother = static_cast<AliMCParticle*>(fMCEvent->GetTrack(imother));
1755         if(McMother->PdgCode() == 111){
1756                 fHistoTrueConvPi0GammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1757         }
1758
1759 }
1760 //________________________________________________________________________
1761 void AliAnalysisTaskGammaConvCalo::ProcessAODMCParticles()
1762 {
1763
1764         TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1765
1766         // Loop over all primary MC particle
1767         for(Int_t i = 0; i < AODMCTrackArray->GetEntriesFast(); i++) {
1768
1769                 AliAODMCParticle* particle = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(i));
1770                 if (!particle) continue;
1771                 if (!particle->IsPrimary()) continue;
1772
1773                 Int_t isMCFromMBHeader = -1;
1774                 if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0){
1775                         isMCFromMBHeader = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent);
1776                         if(isMCFromMBHeader == 0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1777                 }
1778
1779                 if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(particle->Phi(),fEventPlaneAngle,kFALSE)) continue;
1780                 if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(particle,AODMCTrackArray,kFALSE)){
1781                         fHistoMCAllGammaPt[fiCut]->Fill(particle->Pt()); // All MC Gamma
1782                                 if(particle->GetMother() >-1){ // Meson Decay Gamma
1783                                         switch((static_cast<AliAODMCParticle*>(AODMCTrackArray->At(particle->GetMother())))->GetPdgCode()){
1784                                         case 111: // Pi0
1785                                                 fHistoMCDecayGammaPi0Pt[fiCut]->Fill(particle->Pt());
1786                                                 break;
1787                                         case 113: // Rho0
1788                                                 fHistoMCDecayGammaRhoPt[fiCut]->Fill(particle->Pt());
1789                                                 break;
1790                                         case 221: // Eta
1791                                                 fHistoMCDecayGammaEtaPt[fiCut]->Fill(particle->Pt());
1792                                                 break;
1793                                         case 223: // Omega
1794                                                 fHistoMCDecayGammaOmegaPt[fiCut]->Fill(particle->Pt());
1795                                                 break;
1796                                         case 331: // Eta'
1797                                                 fHistoMCDecayGammaEtapPt[fiCut]->Fill(particle->Pt());
1798                                                 break;
1799                                         case 333: // Phi
1800                                                 fHistoMCDecayGammaPhiPt[fiCut]->Fill(particle->Pt());
1801                                                 break;
1802                                         case 3212: // Sigma
1803                                                 fHistoMCDecayGammaSigmaPt[fiCut]->Fill(particle->Pt());
1804                                                 break;
1805                                         }
1806                                 }
1807                         }
1808                         if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(particle,AODMCTrackArray,kTRUE)){
1809                                 Double_t rConv = 0;
1810                                 for(Int_t daughterIndex=particle->GetDaughter(0);daughterIndex<=particle->GetDaughter(1);daughterIndex++){
1811                                         AliAODMCParticle *tmpDaughter = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(daughterIndex));
1812                                         if(!tmpDaughter) continue;
1813                                         if(abs(tmpDaughter->GetPdgCode()) == 11){
1814                                                 rConv = sqrt( (tmpDaughter->Xv()*tmpDaughter->Xv()) + (tmpDaughter->Yv()*tmpDaughter->Yv()) );
1815                                         }
1816                                 }
1817                                 fHistoMCConvGammaPt[fiCut]->Fill(particle->Pt());
1818                                 if (fDoPhotonQA > 0){
1819                                         fHistoMCConvGammaR[fiCut]->Fill(rConv);
1820                                         fHistoMCConvGammaEta[fiCut]->Fill(particle->Eta());
1821                                 }
1822                         }
1823                         // Converted MC Gamma
1824                         if(fDoMesonAnalysis){
1825                         if(particle->GetPdgCode() == 310 && fDoMesonQA > 0){
1826                                 Double_t mesonY = 10.;
1827                                 if(particle->E() - particle->Pz() == 0 || particle->E() + particle->Pz() == 0){
1828                                         mesonY=10.-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift();
1829                                 } else {
1830                                         mesonY = 0.5*(TMath::Log((particle->E()+particle->Pz()) / (particle->E()-particle->Pz())))-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift();
1831                                 }
1832                                 Float_t weightedK0s= 1;
1833                                 if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){
1834                                         if (particle->Pt()>0.005){
1835                                                 weightedK0s= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, 0x0, fInputEvent);
1836                                                 //cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
1837                                         }
1838                                 }
1839                                 fHistoMCK0sPt[fiCut]->Fill(particle->Pt(),weightedK0s);
1840                                 fHistoMCK0sWOWeightPt[fiCut]->Fill(particle->Pt());
1841                                 fHistoMCK0sPtY[fiCut]->Fill(particle->Pt(),mesonY,weightedK0s);
1842                         }
1843                         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
1844                                 ->MesonIsSelectedAODMC(particle,AODMCTrackArray,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift())){
1845                                 AliAODMCParticle* daughter0 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(particle->GetDaughter(0)));
1846                                 AliAODMCParticle* daughter1 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(particle->GetDaughter(1)));
1847                                 Float_t weighted= 1;
1848                                 if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){
1849                                         if (particle->Pt()>0.005){
1850                                                 weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, 0x0, fInputEvent);
1851                                                 //                   if(particle->GetPdgCode() == 221){
1852                                                 //                      cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
1853                                                 //                   }
1854                                         }
1855                                 }
1856                                 Double_t mesonY = 10.;
1857                                 if(particle->E() - particle->Pz() == 0 || particle->E() + particle->Pz() == 0){
1858                                         mesonY=10.-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift();
1859                                 } else{
1860                                         mesonY = 0.5*(TMath::Log((particle->E()+particle->Pz()) / (particle->E()-particle->Pz())))-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift();
1861                                 }
1862
1863                                 if(particle->GetPdgCode() == 111){
1864                                         fHistoMCPi0Pt[fiCut]->Fill(particle->Pt(),weighted); // All MC Pi0
1865                                         fHistoMCPi0WOWeightPt[fiCut]->Fill(particle->Pt());
1866                                         if (fDoMesonQA > 0) fHistoMCPi0PtY[fiCut]->Fill(particle->Pt(),mesonY,weighted); // All MC Pi0
1867                                 } else if(particle->GetPdgCode() == 221){
1868                                         fHistoMCEtaPt[fiCut]->Fill(particle->Pt(),weighted); // All MC Eta
1869                                         fHistoMCEtaWOWeightPt[fiCut]->Fill(particle->Pt());
1870                                         if (fDoMesonQA > 0) fHistoMCEtaPtY[fiCut]->Fill(particle->Pt(),mesonY,weighted); // All MC Pi0
1871                                 }
1872
1873                                 // Check the acceptance for both gammas
1874                                 if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(daughter0,AODMCTrackArray,kFALSE) &&
1875                                 ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(daughter1,AODMCTrackArray,kFALSE)  &&
1876                                 ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter0->Phi(),fEventPlaneAngle,kFALSE) &&
1877                                 ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter1->Phi(),fEventPlaneAngle,kFALSE)){
1878
1879                                         if(particle->GetPdgCode() == 111){
1880                                                 fHistoMCPi0InAccPt[fiCut]->Fill(particle->Pt(),weighted); // MC Pi0 with gamma in acc
1881                                         } else if(particle->GetPdgCode() == 221){
1882                                                 fHistoMCEtaInAccPt[fiCut]->Fill(particle->Pt(),weighted); // MC Eta with gamma in acc
1883                                         }
1884                                 }
1885                         }
1886                 }
1887         }
1888
1889 }
1890 //________________________________________________________________________
1891 void AliAnalysisTaskGammaConvCalo::ProcessMCParticles()
1892 {
1893         // Loop over all primary MC particle
1894         for(Int_t i = 0; i < fMCStack->GetNprimary(); i++) {
1895                 TParticle* particle = (TParticle *)fMCStack->Particle(i);
1896                 if (!particle) continue;
1897
1898                 Int_t isMCFromMBHeader = -1;
1899                 if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0){
1900                         isMCFromMBHeader = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent);
1901                         if(isMCFromMBHeader == 0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1902                 }
1903
1904                 if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(particle->Phi(),fEventPlaneAngle,kFALSE)) continue;
1905                 if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(particle,fMCStack,kFALSE)){
1906                         fHistoMCAllGammaPt[fiCut]->Fill(particle->Pt()); // All MC Gamma
1907                         if(particle->GetMother(0) >-1){ // Meson Decay Gamma
1908                                 switch(fMCStack->Particle(particle->GetMother(0))->GetPdgCode()){
1909                                 case 111: // Pi0
1910                                         fHistoMCDecayGammaPi0Pt[fiCut]->Fill(particle->Pt());
1911                                         break;
1912                                 case 113: // Rho0
1913                                         fHistoMCDecayGammaRhoPt[fiCut]->Fill(particle->Pt());
1914                                         break;
1915                                 case 221: // Eta
1916                                         fHistoMCDecayGammaEtaPt[fiCut]->Fill(particle->Pt());
1917                                         break;
1918                                 case 223: // Omega
1919                                         fHistoMCDecayGammaOmegaPt[fiCut]->Fill(particle->Pt());
1920                                         break;
1921                                 case 331: // Eta'
1922                                         fHistoMCDecayGammaEtapPt[fiCut]->Fill(particle->Pt());
1923                                         break;
1924                                 case 333: // Phi
1925                                         fHistoMCDecayGammaPhiPt[fiCut]->Fill(particle->Pt());
1926                                         break;
1927                                 case 3212: // Sigma
1928                                         fHistoMCDecayGammaSigmaPt[fiCut]->Fill(particle->Pt());
1929                                         break;
1930                                 }
1931                         }
1932                 }
1933                 if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(particle,fMCStack,kTRUE)){
1934                         fHistoMCConvGammaPt[fiCut]->Fill(particle->Pt());
1935                         if (fDoPhotonQA > 0){
1936                                 fHistoMCConvGammaR[fiCut]->Fill(((TParticle*)fMCStack->Particle(particle->GetFirstDaughter()))->R());
1937                                 fHistoMCConvGammaEta[fiCut]->Fill(particle->Eta());
1938                         }
1939                 } // Converted MC Gamma
1940                 if(fDoMesonAnalysis){
1941                         if(particle->GetPdgCode() == 310 && fDoMesonQA > 0){
1942                                 Double_t mesonY = 10.;
1943                                 if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){
1944                                         mesonY=10.-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift();
1945                                 } else{
1946                                         mesonY = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())))-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift();
1947                                 }
1948                                 Float_t weightedK0s= 1;
1949                                 if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){
1950                                         if (particle->Pt()>0.005){
1951                                                 weightedK0s= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, fMCStack, fInputEvent);
1952                                                 //cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
1953                                         }
1954                                 }
1955                                 if (fMCStack->IsPhysicalPrimary(i)){
1956                                         fHistoMCK0sPt[fiCut]->Fill(particle->Pt(),weightedK0s);
1957                                         fHistoMCK0sWOWeightPt[fiCut]->Fill(particle->Pt());
1958                                         fHistoMCK0sPtY[fiCut]->Fill(particle->Pt(),mesonY,weightedK0s);
1959                                 }
1960                         }
1961                         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
1962                                 ->MesonIsSelectedMC(particle,fMCStack,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift())){
1963                                 TParticle* daughter0 = (TParticle*)fMCStack->Particle(particle->GetFirstDaughter());
1964                                 TParticle* daughter1 = (TParticle*)fMCStack->Particle(particle->GetLastDaughter());
1965
1966                                 Float_t weighted= 1;
1967                                 if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){
1968                                         if (particle->Pt()>0.005){
1969                                                 weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, fMCStack, fInputEvent);
1970                                                 //                   if(particle->GetPdgCode() == 221){
1971                                                 //                      cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
1972                                                 //                   }
1973                                         }
1974                                 }
1975                                 Double_t mesonY = 10.;
1976                                 if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){
1977                                         mesonY=10.-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift();
1978                                 } else{
1979                                         mesonY = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())))-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift();
1980                                 }
1981
1982                                 if(particle->GetPdgCode() == 111){
1983                                         fHistoMCPi0Pt[fiCut]->Fill(particle->Pt(),weighted); // All MC Pi0
1984                                         fHistoMCPi0WOWeightPt[fiCut]->Fill(particle->Pt());
1985                                         if (fDoMesonQA > 0) fHistoMCPi0PtY[fiCut]->Fill(particle->Pt(),mesonY,weighted); // All MC Pi0
1986                                 } else if(particle->GetPdgCode() == 221){
1987                                         fHistoMCEtaPt[fiCut]->Fill(particle->Pt(),weighted); // All MC Eta
1988                                         fHistoMCEtaWOWeightPt[fiCut]->Fill(particle->Pt());
1989                                         if (fDoMesonQA > 0) fHistoMCEtaPtY[fiCut]->Fill(particle->Pt(),mesonY,weighted); // All MC Pi0
1990                                 }
1991
1992                                 // Check the acceptance for both gammas
1993                                 if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(daughter0,fMCStack,kFALSE) &&
1994                                 ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(daughter1,fMCStack,kFALSE)  &&
1995                                 ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter0->Phi(),fEventPlaneAngle,kFALSE) &&
1996                                 ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter1->Phi(),fEventPlaneAngle,kFALSE)){
1997
1998                                         if(particle->GetPdgCode() == 111){
1999                                                 fHistoMCPi0InAccPt[fiCut]->Fill(particle->Pt(),weighted); // MC Pi0 with gamma in acc
2000                                         } else if(particle->GetPdgCode() == 221){
2001                                                 fHistoMCEtaInAccPt[fiCut]->Fill(particle->Pt(),weighted); // MC Eta with gamma in acc
2002                                         }
2003                                 }
2004                         }
2005                 }
2006         }
2007
2008         if (fDoMesonQA){
2009                 for(Int_t i = fMCStack->GetNprimary(); i < fMCStack->GetNtrack(); i++) {
2010                         TParticle* particle = (TParticle *)fMCStack->Particle(i);
2011                         if (!particle) continue;
2012
2013                         Int_t isMCFromMBHeader = -1;
2014                         if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0){
2015                                 isMCFromMBHeader = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent);
2016                                 if(isMCFromMBHeader == 0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
2017                         }
2018
2019                         if(fDoMesonAnalysis){
2020                                 if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelectedMC(particle,fMCStack,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift())){
2021                                         Float_t weighted= 1;
2022                                         if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){
2023                                                 if (particle->Pt()>0.005){
2024                                                         weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, fMCStack, fInputEvent);
2025               //                   if(particle->GetPdgCode() == 221){
2026               //                      cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
2027               //                   }
2028                                                 }
2029                                         }
2030
2031                                         if(particle->GetPdgCode() == 111){
2032                                                 Int_t pdgCode = ((TParticle*)fMCStack->Particle( particle->GetFirstMother() ))->GetPdgCode();
2033                                                 Int_t source = GetSourceClassification(111,pdgCode);
2034                                                 fHistoMCSecPi0PtvsSource[fiCut]->Fill(particle->Pt(),source,weighted); // All MC Pi0
2035                                                 fHistoMCSecPi0Source[fiCut]->Fill(pdgCode);
2036                                         } else if(particle->GetPdgCode() == 221){
2037                                                 Int_t pdgCode = ((TParticle*)fMCStack->Particle( particle->GetFirstMother() ))->GetPdgCode();
2038                                                 fHistoMCSecEtaPt[fiCut]->Fill(particle->Pt(),weighted); // All MC Pi0
2039                                                 fHistoMCSecEtaSource[fiCut]->Fill(pdgCode);
2040                                         }
2041                                 }
2042                         }
2043                 }
2044         }
2045 }
2046
2047 //________________________________________________________________________
2048 void AliAnalysisTaskGammaConvCalo::PhotonTagging(){
2049
2050         // Conversion Gammas
2051         if(fGammaCandidates->GetEntries()>0){
2052
2053                 for(Int_t firstGammaIndex=0;firstGammaIndex<fGammaCandidates->GetEntries()-1;firstGammaIndex++){
2054
2055                         // get conversion photon
2056                         AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(firstGammaIndex));
2057                         if (gamma0==NULL) continue;
2058
2059                         TLorentzVector photonVector;
2060                         photonVector.SetPxPyPzE(gamma0->GetPx(),gamma0->GetPy(),gamma0->GetPz(),gamma0->GetPhotonP());
2061
2062                         Bool_t btagpi0 = 0;
2063                         Bool_t btageta = 0;
2064
2065                         // loop over clusters
2066                         for(Int_t secondGammaIndex = 0; secondGammaIndex<fClusterCandidates->GetEntries(); ++secondGammaIndex) {
2067
2068                                 AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fClusterCandidates->At(secondGammaIndex));
2069                                 if (gamma1==NULL) continue;
2070
2071                                 TLorentzVector clusterVector;
2072                                 clusterVector.SetPxPyPzE(gamma1->GetPx(),gamma1->GetPy(),gamma1->GetPz(),gamma1->GetPhotonP());
2073
2074                                 // do the tagging
2075                                 TLorentzVector pairVector = photonVector+clusterVector;
2076
2077                                 // see if pi0?
2078                                 if((pairVector.M() > 0.11 && pairVector.M() < 0.15)){
2079                                         btagpi0 = 1;
2080                                 }
2081                                 // or eta
2082                                 if((pairVector.M() > 0.50 && pairVector.M() < 0.6)){
2083                                         btageta = 1;
2084                                 }
2085                         }// end loop over clusters
2086
2087                         if(btagpi0 && btageta)
2088                                 fHistoConvGammaTagged[fiCut]->Fill(photonVector.Pt());
2089                         else if(btagpi0 && !btageta)
2090                                 fHistoConvGammaPi0Tagged[fiCut]->Fill(photonVector.Pt());
2091                         else if(btageta && !btagpi0)
2092                                 fHistoConvGammaEtaTagged[fiCut]->Fill(photonVector.Pt());
2093                         else
2094                                 fHistoConvGammaUntagged[fiCut]->Fill(photonVector.Pt());
2095
2096                 }// end loop over gammas
2097         }// end if
2098         return;
2099 }
2100
2101 //________________________________________________________________________
2102 void AliAnalysisTaskGammaConvCalo::CalculatePi0Candidates(){
2103
2104         // Conversion Gammas
2105         if(fGammaCandidates->GetEntries()>0){
2106
2107                 // vertex
2108                 Double_t vertex[3] = {0};
2109                 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
2110
2111                 for(Int_t firstGammaIndex=0;firstGammaIndex<fGammaCandidates->GetEntries();firstGammaIndex++){
2112                         AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(firstGammaIndex));
2113                         if (gamma0==NULL) continue;
2114
2115                         for(Int_t secondGammaIndex=0;secondGammaIndex<fClusterCandidates->GetEntries();secondGammaIndex++){
2116
2117                                 AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fClusterCandidates->At(secondGammaIndex));
2118                                 if (gamma1==NULL) continue;
2119
2120                                 AliAODConversionMother *pi0cand = new AliAODConversionMother(gamma0,gamma1);
2121                                 pi0cand->SetLabels(firstGammaIndex,secondGammaIndex);
2122
2123                                 if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(pi0cand,kTRUE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift()))){
2124                                         fHistoMotherInvMassPt[fiCut]->Fill(pi0cand->M(),pi0cand->Pt());
2125
2126                                         // fill new histograms
2127                                         fHistoPhotonPairAll[fiCut]->Fill(pi0cand->M(),pi0cand->Pt());
2128                                         fHistoPhotonPairAllGam[fiCut]->Fill(pi0cand->M(),gamma0->Pt());
2129
2130                                         if(pi0cand->GetAlpha()<0.1)
2131                                                 fHistoMotherInvMassEalpha[fiCut]->Fill(pi0cand->M(),pi0cand->E());
2132
2133                                         if (fDoMesonQA > 0){
2134                                                 if ( pi0cand->M() > 0.05 && pi0cand->M() < 0.17){
2135                                                         fHistoMotherPi0PtY[fiCut]->Fill(pi0cand->Pt(),pi0cand->Rapidity()-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift());
2136                                                         fHistoMotherPi0PtAlpha[fiCut]->Fill(pi0cand->Pt(),pi0cand->GetAlpha());
2137                                                         fHistoMotherPi0PtOpenAngle[fiCut]->Fill(pi0cand->Pt(),pi0cand->GetOpeningAngle());
2138                                                 }
2139                                                 if ( pi0cand->M() > 0.45 && pi0cand->M() < 0.65){
2140                                                         fHistoMotherEtaPtY[fiCut]->Fill(pi0cand->Pt(),pi0cand->Rapidity()-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift());
2141                                                         fHistoMotherEtaPtAlpha[fiCut]->Fill(pi0cand->Pt(),pi0cand->GetAlpha());
2142                                                         fHistoMotherEtaPtOpenAngle[fiCut]->Fill(pi0cand->Pt(),pi0cand->GetOpeningAngle());
2143                                                 }
2144                                         }
2145                                         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->DoBGCalculation()){
2146                                                 Int_t zbin = 0;
2147                                                 Int_t mbin = 0;
2148
2149                                                 if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->BackgroundHandlerType() == 0){
2150                                                         zbin = fBGHandler[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
2151                                                         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2152                                                                 mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
2153                                                         } else {
2154                                                                 mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries());
2155                                                         }
2156                                                 } else{
2157                                                         zbin = fBGHandlerRP[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
2158                                                         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2159                                                                 mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
2160                                                         } else {
2161                                                                 mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries());
2162                                                         }
2163                                                 }
2164                                                 Double_t sparesFill[4] = {pi0cand->M(),pi0cand->Pt(),(Double_t)zbin,(Double_t)mbin};
2165                                                 fSparseMotherInvMassPtZM[fiCut]->Fill(sparesFill,1);
2166                                         }
2167
2168
2169                                         if(fIsMC){
2170                                                 if(fInputEvent->IsA()==AliESDEvent::Class())
2171                                                         ProcessTrueMesonCandidates(pi0cand,gamma0,gamma1);
2172                                                 if(fInputEvent->IsA()==AliAODEvent::Class())
2173                                                         ProcessTrueMesonCandidatesAOD(pi0cand,gamma0,gamma1);
2174                                         }
2175                                         if (fDoMesonQA == 2){
2176                                                 fInvMass = pi0cand->M();
2177                                                 fPt  = pi0cand->Pt();
2178                                                 if (abs(gamma0->GetDCAzToPrimVtx()) < abs(gamma1->GetDCAzToPrimVtx())){
2179                                                         fDCAzGammaMin = gamma0->GetDCAzToPrimVtx();
2180                                                         fDCAzGammaMax = gamma1->GetDCAzToPrimVtx();
2181                                                 } else {
2182                                                         fDCAzGammaMin = gamma1->GetDCAzToPrimVtx();
2183                                                         fDCAzGammaMax = gamma0->GetDCAzToPrimVtx();
2184                                                 }
2185                                                 fCharFlag = pi0cand->GetMesonQuality();
2186                                                 //                   cout << "gamma 0: " << gamma0->GetV0Index()<< "\t" << gamma0->GetPx() << "\t" << gamma0->GetPy() << "\t" <<  gamma0->GetPz() << "\t" << endl;
2187                                                 //                   cout << "gamma 1: " << gamma1->GetV0Index()<< "\t"<< gamma1->GetPx() << "\t" << gamma1->GetPy() << "\t" <<  gamma1->GetPz() << "\t" << endl;
2188                                                 //                    cout << "pi0: "<<fInvMass << "\t" << fPt <<"\t" << fDCAzGammaMin << "\t" << fDCAzGammaMax << "\t" << (Int_t)fCharFlag << "\t" << (Int_t)fCharMesonMCInfo <<endl;
2189                                                 if (fIsHeavyIon == 1 && fPt > 0.399 && fPt < 20. ) {
2190                                                 if (fInvMass > 0.08 && fInvMass < 0.2) fTreeMesonsInvMassPtDcazMinDcazMaxFlag[fiCut]->Fill();
2191                                                 if ((fInvMass > 0.45 && fInvMass < 0.6) &&  (fPt > 0.999 && fPt < 20.) )fTreeMesonsInvMassPtDcazMinDcazMaxFlag[fiCut]->Fill();
2192                                                 } else if (fPt > 0.299 && fPt < 20. )  {
2193                                                 if ( (fInvMass > 0.08 && fInvMass < 0.2) || (fInvMass > 0.45 && fInvMass < 0.6)) fTreeMesonsInvMassPtDcazMinDcazMaxFlag[fiCut]->Fill();
2194                                                 }
2195                                         }
2196                                 }
2197                                 delete pi0cand;
2198                                 pi0cand=0x0;
2199                         }
2200                 }
2201         }
2202 }
2203 //______________________________________________________________________
2204 void AliAnalysisTaskGammaConvCalo::ProcessTrueMesonCandidates(AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1)
2205 {
2206         // Process True Mesons
2207         AliStack *MCStack = fMCEvent->Stack();
2208         fCharMesonMCInfo = 0;
2209         if(TrueGammaCandidate0->GetV0Index()<fInputEvent->GetNumberOfV0s()){
2210                 Bool_t isTruePi0 = kFALSE;
2211                 Bool_t isTrueEta = kFALSE;
2212                 Bool_t isTruePi0Dalitz = kFALSE;
2213                 Bool_t isTrueEtaDalitz = kFALSE;
2214                 Bool_t gamma0DalitzCand = kFALSE;
2215                 Bool_t gamma1DalitzCand = kFALSE;
2216                 Int_t gamma0MCLabel = TrueGammaCandidate0->GetMCParticleLabel(MCStack);
2217                 Int_t gamma0MotherLabel = -1;
2218                 if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
2219                 // Daughters Gamma 0
2220                 TParticle * negativeMC = (TParticle*)TrueGammaCandidate0->GetNegativeMCDaughter(MCStack);
2221                 TParticle * positiveMC = (TParticle*)TrueGammaCandidate0->GetPositiveMCDaughter(MCStack);
2222                 TParticle * gammaMC0 = (TParticle*)MCStack->Particle(gamma0MCLabel);
2223                 if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
2224                         if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
2225                                 if(gammaMC0->GetPdgCode() == 22){ // ... with Gamma Mother
2226                                         gamma0MotherLabel=gammaMC0->GetFirstMother();
2227                                 }
2228                         }
2229                         if(gammaMC0->GetPdgCode() ==111){ // Dalitz candidate
2230                                 gamma0DalitzCand = kTRUE;
2231                                 gamma0MotherLabel=-111;
2232                         }
2233                         if(gammaMC0->GetPdgCode() ==221){ // Dalitz candidate
2234                                 gamma0DalitzCand = kTRUE;
2235                                 gamma0MotherLabel=-221;
2236                         }
2237                 }
2238                 }
2239                 if(TrueGammaCandidate1->GetV0Index()<fInputEvent->GetNumberOfV0s()){
2240                         Int_t gamma1MCLabel = TrueGammaCandidate1->GetMCParticleLabel(MCStack);
2241                         Int_t gamma1MotherLabel = -1;
2242                         if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
2243                                 // Daughters Gamma 1
2244                                 TParticle * negativeMC = (TParticle*)TrueGammaCandidate1->GetNegativeMCDaughter(MCStack);
2245                                 TParticle * positiveMC = (TParticle*)TrueGammaCandidate1->GetPositiveMCDaughter(MCStack);
2246                                 TParticle * gammaMC1 = (TParticle*)MCStack->Particle(gamma1MCLabel);
2247                                 if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
2248                                         if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
2249                                                 if(gammaMC1->GetPdgCode() == 22){ // ... with Gamma Mother
2250                                                 gamma1MotherLabel=gammaMC1->GetFirstMother();
2251                                                 }
2252                                         }
2253                                         if(gammaMC1->GetPdgCode() ==111 ){ // Dalitz candidate
2254                                                 gamma1DalitzCand = kTRUE;
2255                                                 gamma1MotherLabel=-111;
2256                                         }
2257                                         if(gammaMC1->GetPdgCode() ==221){ // Dalitz candidate
2258                                                 gamma1DalitzCand = kTRUE;
2259                                                 gamma1MotherLabel=-221;
2260                                         }
2261                                 }
2262                         }
2263                         if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
2264                                 if(((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetPdgCode() == 111){
2265                                         isTruePi0=kTRUE;
2266                                 }
2267                                 if(((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetPdgCode() == 221){
2268                                         isTrueEta=kTRUE;
2269                                 }
2270                         }
2271
2272                         //Identify Dalitz candidate
2273                         if (gamma1DalitzCand || gamma0DalitzCand){
2274                                 if (gamma0DalitzCand && gamma0MCLabel >=0 && gamma0MCLabel==gamma1MotherLabel){
2275                                         if (gamma0MotherLabel == -111) isTruePi0Dalitz = kTRUE;
2276                                         if (gamma0MotherLabel == -221) isTrueEtaDalitz = kTRUE;
2277                                 }
2278                                 if (gamma1DalitzCand && gamma1MCLabel >=0 && gamma1MCLabel==gamma0MotherLabel){
2279                                         if (gamma1MotherLabel == -111) isTruePi0Dalitz = kTRUE;
2280                                         if (gamma1MotherLabel == -221) isTrueEtaDalitz = kTRUE;
2281                                 }
2282                         }
2283
2284
2285                         if(isTruePi0 || isTrueEta){// True Pion or Eta
2286                                 fHistoTrueMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2287                                 if (fDoMesonQA > 0){
2288                                         if (isTruePi0){
2289                                                 if ( Pi0Candidate->M() > 0.05 && Pi0Candidate->M() < 0.17){
2290                                                         fHistoTruePi0PtY[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->Rapidity()-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift());
2291                                                         fHistoTruePi0PtAlpha[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetAlpha());
2292                                                         fHistoTruePi0PtOpenAngle[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetOpeningAngle());
2293                                                 }
2294                                         } else if (isTrueEta){
2295                                                 if ( Pi0Candidate->M() > 0.45 && Pi0Candidate->M() < 0.65){
2296                                                         fHistoTrueEtaPtY[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->Rapidity()-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift());
2297                                                         fHistoTrueEtaPtAlpha[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetAlpha());
2298                                                         fHistoTrueEtaPtOpenAngle[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetOpeningAngle());
2299                                                 }
2300                                         }
2301                                 }
2302                                 if(gamma0MotherLabel >= MCStack->GetNprimary()){ // Secondary Meson
2303                                         Int_t secMotherLabel = ((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetMother(0);
2304                                         Float_t weightedSec= 1;
2305                                         if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(secMotherLabel, fMCStack, fInputEvent) && MCStack->Particle(secMotherLabel)->GetPdgCode()==310){
2306                                                 weightedSec= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),secMotherLabel, fMCStack, fInputEvent)/2.; //invariant mass is additive thus the weight for the daughters has to be devide by two for the K0s at a certain pt
2307                                                 //cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
2308                                         }
2309                                         fHistoTrueSecondaryMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2310                                         fCharMesonMCInfo = 2;
2311                                         if (secMotherLabel >-1){
2312                                                 if(MCStack->Particle(secMotherLabel)->GetPdgCode()==310){
2313                                                         fCharMesonMCInfo = 4;
2314                                                         fHistoTrueSecondaryMotherFromK0sInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2315                                                         if (fDoMesonQA > 0)fHistoTrueK0sWithPi0DaughterMCPt[fiCut]->Fill(MCStack->Particle(secMotherLabel)->Pt());
2316                                                 }
2317                                                 if(MCStack->Particle(secMotherLabel)->GetPdgCode()==221){
2318                                                 fCharMesonMCInfo = 3;
2319                                                 fHistoTrueSecondaryMotherFromEtaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2320                                                 if (fDoMesonQA > 0)fHistoTrueEtaWithPi0DaughterMCPt[fiCut]->Fill(MCStack->Particle(secMotherLabel)->Pt());
2321                                                 }
2322                                                 if(MCStack->Particle(secMotherLabel)->GetPdgCode()==3122){
2323                                                         fCharMesonMCInfo = 7;
2324                                                         fHistoTrueSecondaryMotherFromLambdaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2325                                                         if (fDoMesonQA > 0)fHistoTrueLambdaWithPi0DaughterMCPt[fiCut]->Fill(MCStack->Particle(secMotherLabel)->Pt());
2326                                                 }
2327                                         }
2328                                 } else { // Only primary pi0 for efficiency calculation
2329                                         fCharMesonMCInfo = 6;
2330                                         Float_t weighted= 1;
2331                                         if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(gamma1MotherLabel, fMCStack, fInputEvent)){
2332                                                 if (((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt()>0.005){
2333                                                         weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),gamma1MotherLabel, fMCStack, fInputEvent);
2334                                                         //                      cout << "rec \t " <<gamma1MotherLabel << "\t" <<  weighted << endl;
2335                                                 }
2336                                         }
2337                                         fHistoTruePrimaryMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weighted);
2338                                         fHistoTruePrimaryMotherW0WeightingInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2339                                         fProfileTruePrimaryMotherWeightsInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weighted);
2340
2341
2342                                         if (fDoMesonQA > 0){
2343                                                 if(isTruePi0){ // Only primary pi0 for resolution
2344                                                         fHistoTruePrimaryPi0MCPtResolPt[fiCut]->Fill(((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt(),(Pi0Candidate->Pt()-((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt())/((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt(),weighted);
2345                                                 }
2346                                                 if (isTrueEta){ // Only primary eta for resolution
2347                                                         fHistoTruePrimaryEtaMCPtResolPt[fiCut]->Fill(((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt(),(Pi0Candidate->Pt()-((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt())/((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt(),weighted);
2348                                                 }
2349                                         }
2350                                 }
2351                         } else if(!isTruePi0 && !isTrueEta){ // Background
2352                                 if (fDoMesonQA > 0){
2353                                         if(gamma0MotherLabel>-1 && gamma1MotherLabel>-1){ // Both Tracks are Photons and have a mother but not Pi0 or Eta
2354                                                 fHistoTrueBckGGInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2355                                                 fCharMesonMCInfo = 1;
2356                                         } else { // No photon or without mother
2357                                                 fHistoTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2358                                         }
2359                                 }
2360                                 if( isTruePi0Dalitz || isTrueEtaDalitz ){
2361                                         // Dalitz
2362                                         fCharMesonMCInfo = 5;
2363                                         fHistoTrueMotherDalitzInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2364                                 } else if (gamma0DalitzCand || gamma1DalitzCand){
2365                                         if (fDoMesonQA > 0)fHistoTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2366                                 }
2367                         }
2368                 }
2369         }
2370 }
2371 //______________________________________________________________________
2372 void AliAnalysisTaskGammaConvCalo::ProcessTrueMesonCandidatesAOD(AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1)
2373 {
2374
2375         // Process True Mesons
2376         TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
2377         Bool_t isTruePi0 = kFALSE;
2378         Bool_t isTrueEta = kFALSE;
2379         Bool_t isTruePi0Dalitz = kFALSE;
2380         Bool_t isTrueEtaDalitz = kFALSE;
2381         Bool_t gamma0DalitzCand = kFALSE;
2382         Bool_t gamma1DalitzCand = kFALSE;
2383
2384         AliAODMCParticle *positiveMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate0->GetMCLabelPositive()));
2385         AliAODMCParticle *negativeMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate0->GetMCLabelNegative()));
2386
2387         fCharMesonMCInfo = 0;
2388         Int_t gamma0MCLabel = -1;
2389         Int_t gamma0MotherLabel = -1;
2390         if(!positiveMC||!negativeMC)
2391                 return;
2392
2393         if(positiveMC->GetMother()>-1&&(negativeMC->GetMother() == positiveMC->GetMother())){
2394                 gamma0MCLabel = positiveMC->GetMother();
2395         }
2396
2397         if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
2398                 // Daughters Gamma 0
2399                 AliAODMCParticle * gammaMC0 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma0MCLabel));
2400                 if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
2401                         if(((positiveMC->GetMCProcessCode())) == 5 && ((negativeMC->GetMCProcessCode())) == 5){ // ... From Conversion ...
2402                                 if(gammaMC0->GetPdgCode() == 22){ // ... with Gamma Mother
2403                                 gamma0MotherLabel=gammaMC0->GetMother();
2404                                 }
2405                         }
2406                         if(gammaMC0->GetPdgCode() ==111){ // Dalitz candidate
2407                                 gamma0DalitzCand = kTRUE;
2408                                 gamma0MotherLabel=-111;
2409                         }
2410                         if(gammaMC0->GetPdgCode() ==221){ // Dalitz candidate
2411                                 gamma0DalitzCand = kTRUE;
2412                                 gamma0MotherLabel=-221;
2413                         }
2414                 }
2415         }
2416         positiveMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate1->GetMCLabelPositive()));
2417         negativeMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate1->GetMCLabelNegative()));
2418
2419         Int_t gamma1MCLabel = -1;
2420         Int_t gamma1MotherLabel = -1;
2421         if(!positiveMC||!negativeMC)
2422                 return;
2423
2424         if(positiveMC->GetMother()>-1&&(negativeMC->GetMother() == positiveMC->GetMother())){
2425                 gamma1MCLabel = positiveMC->GetMother();
2426         }
2427         if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
2428                 // Daughters Gamma 1
2429                 AliAODMCParticle * gammaMC1 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MCLabel));
2430                 if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
2431                         if(((positiveMC->GetMCProcessCode())) == 5 && ((negativeMC->GetMCProcessCode())) == 5){ // ... From Conversion ...
2432                                 if(gammaMC1->GetPdgCode() == 22){ // ... with Gamma Mother
2433                                 gamma1MotherLabel=gammaMC1->GetMother();
2434                                 }
2435                         }
2436                         if(gammaMC1->GetPdgCode() ==111 ){ // Dalitz candidate
2437                                 gamma1DalitzCand = kTRUE;
2438                                 gamma1MotherLabel=-111;
2439                         }
2440                         if(gammaMC1->GetPdgCode() ==221){ // Dalitz candidate
2441                                 gamma1DalitzCand = kTRUE;
2442                                 gamma1MotherLabel=-221;
2443                         }
2444                 }
2445         }
2446         if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
2447                 if(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->GetPdgCode() == 111){
2448                         isTruePi0=kTRUE;
2449                 }
2450                 if(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->GetPdgCode() == 221){
2451                         isTrueEta=kTRUE;
2452                 }
2453         }
2454
2455         //Identify Dalitz candidate
2456         if (gamma1DalitzCand || gamma0DalitzCand){
2457                 if (gamma0DalitzCand && gamma0MCLabel >=0 && gamma0MCLabel==gamma1MotherLabel){
2458                         if (gamma0MotherLabel == -111) isTruePi0Dalitz = kTRUE;
2459                         if (gamma0MotherLabel == -221) isTrueEtaDalitz = kTRUE;
2460                 }
2461                 if (gamma1DalitzCand && gamma1MCLabel >=0 && gamma1MCLabel==gamma0MotherLabel){
2462                         if (gamma1MotherLabel == -111) isTruePi0Dalitz = kTRUE;
2463                         if (gamma1MotherLabel == -221) isTrueEtaDalitz = kTRUE;
2464                 }
2465         }
2466
2467         if(isTruePi0 || isTrueEta){// True Pion or Eta
2468                 fHistoTrueMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2469                 if (fDoMesonQA > 0){
2470                         if (isTruePi0){
2471                                 if ( Pi0Candidate->M() > 0.05 && Pi0Candidate->M() < 0.17){
2472                                 fHistoTruePi0PtY[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->Rapidity()-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift());
2473                                 fHistoTruePi0PtAlpha[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetAlpha());
2474                                 fHistoTruePi0PtOpenAngle[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetOpeningAngle());
2475                                 }
2476                         } else if (isTrueEta){
2477                                 if ( Pi0Candidate->M() > 0.45 && Pi0Candidate->M() < 0.65){
2478                                 fHistoTrueEtaPtY[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->Rapidity()-((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift());
2479                                 fHistoTrueEtaPtAlpha[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetAlpha());
2480                                 fHistoTrueEtaPtOpenAngle[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetOpeningAngle());
2481                                 }
2482                         }
2483                 }
2484                 if(!(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma0MotherLabel))->IsPrimary())){ // Secondary Meson
2485                         Int_t secMotherLabel = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->GetMother();
2486                         Float_t weightedSec= 1;
2487                         if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(secMotherLabel, 0x0, fInputEvent) && static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->GetPdgCode()==310){
2488                                 weightedSec= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),secMotherLabel, 0x0, fInputEvent)/2.; //invariant mass is additive thus the weight for the daughters has to be devide by two for the K0s at a certain pt
2489                                 //cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
2490                         }
2491                         fHistoTrueSecondaryMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2492                         fCharMesonMCInfo = 2;
2493                         if (secMotherLabel >-1){
2494                                 if(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->GetPdgCode()==310){
2495                                         fCharMesonMCInfo = 4;
2496                                         fHistoTrueSecondaryMotherFromK0sInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2497                                         if (fDoMesonQA > 0)fHistoTrueK0sWithPi0DaughterMCPt[fiCut]->Fill(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->Pt());
2498                                 }
2499                                 if(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->GetPdgCode()==221){
2500                                         fCharMesonMCInfo = 3;
2501                                         fHistoTrueSecondaryMotherFromEtaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2502                                         if (fDoMesonQA > 0)fHistoTrueEtaWithPi0DaughterMCPt[fiCut]->Fill(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->Pt());
2503                                 }
2504                                 if(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->GetPdgCode()==3122){
2505                                         fCharMesonMCInfo = 7;
2506                                         fHistoTrueSecondaryMotherFromLambdaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2507                                         if (fDoMesonQA > 0)fHistoTrueLambdaWithPi0DaughterMCPt[fiCut]->Fill(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->Pt());
2508                                 }
2509                         }
2510                 }else{ // Only primary pi0 for efficiency calculation
2511                         Float_t weighted= 1;
2512                         fCharMesonMCInfo = 6;
2513                         if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(gamma1MotherLabel, 0x0, fInputEvent)){
2514                                 if (static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt()>0.005){
2515                                 weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),gamma1MotherLabel, 0x0, fInputEvent);
2516                                 //                      cout << "rec \t " <<gamma1MotherLabel << "\t" <<  weighted << endl;
2517                                 }
2518                         }
2519                         fHistoTruePrimaryMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weighted);
2520                         fHistoTruePrimaryMotherW0WeightingInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2521                         fProfileTruePrimaryMotherWeightsInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weighted);
2522
2523                         if (fDoMesonQA > 0){
2524                                 if(isTruePi0){ // Only primary pi0 for resolution
2525                                         fHistoTruePrimaryPi0MCPtResolPt[fiCut]->Fill(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),
2526                                                                                                                         (Pi0Candidate->Pt()-static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt())/static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),weighted);
2527
2528                                 }
2529                                 if (isTrueEta){ // Only primary eta for resolution
2530                                         fHistoTruePrimaryEtaMCPtResolPt[fiCut]->Fill(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),
2531                                                                                                                         (Pi0Candidate->Pt()-static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt())/static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),weighted);
2532                                 }
2533                         }
2534                 }
2535         } else if(!isTruePi0 && !isTrueEta) { // Background
2536                 if (fDoMesonQA > 0){
2537                         if(gamma0MotherLabel>-1 && gamma1MotherLabel>-1){ // Both Tracks are Photons and have a mother but not Pi0 or Eta
2538                                 fHistoTrueBckGGInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2539                                 fCharMesonMCInfo = 1;
2540                         } else { // No photon or without mother
2541                                 fHistoTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2542                         }
2543                 }
2544                 if( isTruePi0Dalitz || isTrueEtaDalitz ){
2545                         // Dalitz
2546                         fCharMesonMCInfo = 5;
2547                         fHistoTrueMotherDalitzInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2548                 } else if (gamma0DalitzCand || gamma1DalitzCand){
2549                         if (fDoMesonQA > 0)fHistoTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2550                 }
2551         }
2552 }
2553
2554 //________________________________________________________________________
2555 void AliAnalysisTaskGammaConvCalo::CalculateBackground(){
2556
2557         Int_t zbin= fBGClusHandler[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
2558         Int_t mbin = 0;
2559
2560         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2561                 mbin = fBGClusHandler[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
2562         } else {
2563                 mbin = fBGClusHandler[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries());
2564         }
2565
2566         AliGammaConversionAODBGHandler::GammaConversionVertex *bgEventVertex = NULL;
2567         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2568                 for(Int_t nEventsInBG=0;nEventsInBG<fBGClusHandler[fiCut]->GetNBGEvents();nEventsInBG++){
2569                         AliGammaConversionAODVector *previousEventV0s = fBGClusHandler[fiCut]->GetBGGoodV0s(zbin,mbin,nEventsInBG);
2570                         if(fMoveParticleAccordingToVertex == kTRUE){
2571                                 bgEventVertex = fBGClusHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG);
2572                         }
2573
2574                         for(Int_t iCurrent=0;iCurrent<fGammaCandidates->GetEntries();iCurrent++){
2575                                 AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGammaCandidates->At(iCurrent));
2576                                 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2577                                         AliAODConversionPhoton previousGoodV0 = (AliAODConversionPhoton)(*(previousEventV0s->at(iPrevious)));
2578                                         if(fMoveParticleAccordingToVertex == kTRUE){
2579                                                 MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
2580                                         }
2581                                         if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->GetInPlaneOutOfPlaneCut() != 0){
2582                                                 RotateParticleAccordingToEP(&previousGoodV0,bgEventVertex->fEP,fEventPlaneAngle);
2583                                         }
2584
2585                                         AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&previousGoodV0);
2586                                         backgroundCandidate->CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
2587                                         if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
2588                                                 ->MesonIsSelected(backgroundCandidate,kFALSE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift()))){
2589                                                 fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
2590                                                 Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin};
2591                                                 fSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
2592                                         }
2593                                         delete backgroundCandidate;
2594                                         backgroundCandidate = 0x0;
2595                                 }
2596                         }
2597                 }
2598         } else {
2599                 for(Int_t nEventsInBG=0;nEventsInBG <fBGClusHandler[fiCut]->GetNBGEvents();nEventsInBG++){
2600                         AliGammaConversionAODVector *previousEventV0s = fBGClusHandler[fiCut]->GetBGGoodV0s(zbin,mbin,nEventsInBG);
2601                         if(previousEventV0s){
2602                                 if(fMoveParticleAccordingToVertex == kTRUE){
2603                                         bgEventVertex = fBGClusHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG);
2604                                 }
2605                                 for(Int_t iCurrent=0;iCurrent<fGammaCandidates->GetEntries();iCurrent++){
2606                                         AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGammaCandidates->At(iCurrent));
2607                                         for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2608
2609                                                 AliAODConversionPhoton previousGoodV0 = (AliAODConversionPhoton)(*(previousEventV0s->at(iPrevious)));
2610
2611                                                 if(fMoveParticleAccordingToVertex == kTRUE){
2612                                                         MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
2613                                                 }
2614                                                 if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->GetInPlaneOutOfPlaneCut() != 0){
2615                                                         RotateParticleAccordingToEP(&previousGoodV0,bgEventVertex->fEP,fEventPlaneAngle);
2616                                                 }
2617
2618
2619                                                 AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&previousGoodV0);
2620                                                 backgroundCandidate->CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
2621                                                 if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift()))){
2622                                                         fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
2623                                                         Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin};
2624                                                         fSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
2625                                                 }
2626                                                 delete backgroundCandidate;
2627                                                 backgroundCandidate = 0x0;
2628                                         }
2629                                 }
2630                         }
2631                 }
2632         }
2633 }
2634
2635 //________________________________________________________________________
2636 void AliAnalysisTaskGammaConvCalo::CalculateBackgroundRP(){
2637
2638         Int_t zbin= fBGHandlerRP[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
2639         Int_t mbin = 0;
2640         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2641                 mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
2642         } else {
2643                 mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries());
2644         }
2645
2646
2647         //Rotation Method
2648         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseRotationMethod()){
2649                 // Correct for the number of rotations
2650                 // BG is for rotation the same, except for factor NRotations
2651                 Double_t weight=1./Double_t(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->GetNumberOfBGEvents());
2652
2653                 for(Int_t firstGammaIndex=0;firstGammaIndex<fGammaCandidates->GetEntries();firstGammaIndex++){
2654
2655                         AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(firstGammaIndex));
2656                         if (gamma0==NULL) continue;
2657                         for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fGammaCandidates->GetEntries();secondGammaIndex++){
2658                                 AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(secondGammaIndex));
2659                                 if (gamma1 == NULL) continue;
2660                                 if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelected(gamma1,fInputEvent))continue;
2661                                 for(Int_t nRandom=0;nRandom<((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->GetNumberOfBGEvents();nRandom++){
2662
2663                                         RotateParticle(gamma1);
2664                                         AliAODConversionMother backgroundCandidate(gamma0,gamma1);
2665                                         backgroundCandidate.CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
2666                                         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
2667                                                 ->MesonIsSelected(&backgroundCandidate,kFALSE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift())){
2668                                                 fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate.M(),backgroundCandidate.Pt());
2669                                                 Double_t sparesFill[4] = {backgroundCandidate.M(),backgroundCandidate.Pt(),(Double_t)zbin,(Double_t)mbin};
2670                                                 fSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,weight);
2671                                         }
2672                                 }
2673                         }
2674                 }
2675         } else {
2676                 // Do Event Mixing
2677                 for(Int_t nEventsInBG=0;nEventsInBG <fBGHandlerRP[fiCut]->GetNBGEvents(fGammaCandidates,fInputEvent);nEventsInBG++){
2678
2679                         AliGammaConversionPhotonVector *previousEventGammas = fBGHandlerRP[fiCut]->GetBGGoodGammas(fGammaCandidates,fInputEvent,nEventsInBG);
2680
2681                         if(previousEventGammas){
2682                                 // test weighted background
2683                                 Double_t weight=1.0;
2684                                 // Correct for the number of eventmixing:
2685                                 // N gammas -> (N-1) + (N-2) +(N-3) ...+ (N-(N-1))  using sum formula sum(i)=N*(N-1)/2  -> N*(N-1)/2
2686                                 // real combinations (since you cannot combine a photon with its own)
2687                                 // but BG leads to N_{a}*N_{b} combinations
2688                                 weight*=0.5*(Double_t(fGammaCandidates->GetEntries()-1))/Double_t(previousEventGammas->size());
2689
2690                                 for(Int_t iCurrent=0;iCurrent<fGammaCandidates->GetEntries();iCurrent++){
2691                                         AliAODConversionPhoton *gamma0 = (AliAODConversionPhoton*)(fGammaCandidates->At(iCurrent));
2692                                         for(UInt_t iPrevious=0;iPrevious<previousEventGammas->size();iPrevious++){
2693
2694                                                 AliAODConversionPhoton *gamma1 = (AliAODConversionPhoton*)(previousEventGammas->at(iPrevious));
2695
2696                                                 AliAODConversionMother backgroundCandidate(gamma0,gamma1);
2697                                                 backgroundCandidate.CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
2698                                                 if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
2699                                                 ->MesonIsSelected(&backgroundCandidate,kFALSE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift())){
2700                                                 fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate.M(),backgroundCandidate.Pt());
2701                                                 Double_t sparesFill[4] = {backgroundCandidate.M(),backgroundCandidate.Pt(),(Double_t)zbin,(Double_t)mbin};
2702                                                 fSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,weight);
2703                                                 }
2704                                         }
2705                                 }
2706                         }
2707                 }
2708         }
2709 }
2710 //________________________________________________________________________
2711 void AliAnalysisTaskGammaConvCalo::RotateParticle(AliAODConversionPhoton *gamma){
2712         Int_t fNDegreesPMBackground= ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->NDegreesRotation();
2713         Double_t nRadiansPM = fNDegreesPMBackground*TMath::Pi()/180;
2714         Double_t rotationValue = fRandom.Rndm()*2*nRadiansPM + TMath::Pi()-nRadiansPM;
2715         gamma->RotateZ(rotationValue);
2716 }
2717
2718 //________________________________________________________________________
2719 void AliAnalysisTaskGammaConvCalo::RotateParticleAccordingToEP(AliAODConversionPhoton *gamma, Double_t previousEventEP, Double_t thisEventEP){
2720
2721         previousEventEP=previousEventEP+TMath::Pi();
2722         thisEventEP=thisEventEP+TMath::Pi();
2723         Double_t rotationValue= thisEventEP-previousEventEP;
2724         gamma->RotateZ(rotationValue);
2725 }
2726
2727 //________________________________________________________________________
2728 void AliAnalysisTaskGammaConvCalo::MoveParticleAccordingToVertex(AliAODConversionPhoton* particle,const AliGammaConversionAODBGHandler::GammaConversionVertex *vertex){
2729         //see header file for documentation
2730
2731         Double_t dx = vertex->fX - fInputEvent->GetPrimaryVertex()->GetX();
2732         Double_t dy = vertex->fY - fInputEvent->GetPrimaryVertex()->GetY();
2733         Double_t dz = vertex->fZ - fInputEvent->GetPrimaryVertex()->GetZ();
2734
2735         Double_t movedPlace[3] = {particle->GetConversionX() - dx,particle->GetConversionY() - dy,particle->GetConversionZ() - dz};
2736         particle->SetConversionPoint(movedPlace);
2737 }
2738 //________________________________________________________________________
2739 void AliAnalysisTaskGammaConvCalo::UpdateEventByEventData(){
2740         //see header file for documentation
2741         if(fGammaCandidates->GetEntries() >0 ){
2742                 if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2743                         fBGHandler[fiCut]->AddEvent(fGammaCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fV0Reader->GetNumberOfPrimaryTracks(),fEventPlaneAngle);
2744                         fBGClusHandler[fiCut]->AddEvent(fClusterCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fV0Reader->GetNumberOfPrimaryTracks(),fEventPlaneAngle);
2745                 } else { // means we use #V0s for multiplicity
2746                         fBGHandler[fiCut]->AddEvent(fGammaCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fGammaCandidates->GetEntries(),fEventPlaneAngle);
2747                         fBGClusHandler[fiCut]->AddEvent(fClusterCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fGammaCandidates->GetEntries(),fEventPlaneAngle);
2748                 }
2749         }
2750 }
2751
2752
2753 //________________________________________________________________________
2754 void AliAnalysisTaskGammaConvCalo::FillPhotonCombinatorialBackgroundHist(AliAODConversionPhoton *TruePhotonCandidate, Int_t pdgCode[])
2755 {
2756         // Combinatorial Bck = 0 ee, 1 ep,i 2 ek, 3 ep, 4 emu, 5 pipi, 6 pik, 7 pip, 8 pimu, 9 kk, 10 kp, 11 kmu, 12 pp, 13 pmu, 14 mumu, 15 Rest
2757         if(pdgCode[0]==11   && pdgCode[1]==11){if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),0);}
2758         else if( (pdgCode[0]==11   && pdgCode[1]==211) || (pdgCode[0]==211  && pdgCode[1]==11) )
2759         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),1);}
2760         else if( (pdgCode[0]==11   && pdgCode[1]==321) || (pdgCode[0]==321  && pdgCode[1]==11) )
2761         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),2);}
2762         else if( (pdgCode[0]==11   && pdgCode[1]==2212) || (pdgCode[0]==2212 && pdgCode[1]==11) )
2763         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),3);}
2764         else if( (pdgCode[0]==11   && pdgCode[1]==13) || (pdgCode[0]==13   && pdgCode[1]==11) )
2765         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),4);}
2766         else if(  pdgCode[0]==211  && pdgCode[1]==211 ){if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),5);}
2767         else if( (pdgCode[0]==211  && pdgCode[1]==321) || (pdgCode[0]==321  && pdgCode[1]==211) )
2768         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),6);}
2769         else if( (pdgCode[0]==211  && pdgCode[1]==2212) || (pdgCode[0]==2212 && pdgCode[1]==211) )
2770         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),7);}
2771         else if( (pdgCode[0]==211  && pdgCode[1]==13) || (pdgCode[0]==13   && pdgCode[1]==211) )
2772         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),8);}
2773         else if(  pdgCode[0]==321  && pdgCode[1]==321 ){if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),9);}
2774         else if( (pdgCode[0]==321  && pdgCode[1]==2212) || (pdgCode[0]==2212 && pdgCode[1]==321) )
2775         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),10);}
2776         else if( (pdgCode[0]==321  && pdgCode[1]==13) || (pdgCode[0]==13   && pdgCode[1]==321) )
2777         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),11);}
2778         else if(  pdgCode[0]==2212   && pdgCode[1]==2212  ){if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),12);}
2779         else if( (pdgCode[0]==2212  && pdgCode[1]==13) || (pdgCode[0]==13   && pdgCode[1]==2212) )
2780         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),13);}
2781         else if(  pdgCode[0]==13   && pdgCode[1]==13  ){if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),14);}
2782         else {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),15);}
2783 }
2784 //________________________________________________________________________
2785 void AliAnalysisTaskGammaConvCalo::RelabelAODPhotonCandidates(Bool_t mode){
2786
2787         // Relabeling For AOD Event
2788         // ESDiD -> AODiD
2789         // MCLabel -> AODMCLabel
2790
2791         if(mode){
2792                 fMCStackPos = new Int_t[fReaderGammas->GetEntries()];
2793                 fMCStackNeg = new Int_t[fReaderGammas->GetEntries()];
2794                 fESDArrayPos = new Int_t[fReaderGammas->GetEntries()];
2795                 fESDArrayNeg = new Int_t[fReaderGammas->GetEntries()];
2796         }
2797
2798         for(Int_t iGamma = 0;iGamma<fReaderGammas->GetEntries();iGamma++){
2799                 AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) fReaderGammas->At(iGamma);
2800                 if(!PhotonCandidate) continue;
2801                 if(!mode){// Back to ESD Labels
2802                 PhotonCandidate->SetMCLabelPositive(fMCStackPos[iGamma]);
2803                 PhotonCandidate->SetMCLabelNegative(fMCStackNeg[iGamma]);
2804                 PhotonCandidate->SetLabelPositive(fESDArrayPos[iGamma]);
2805                 PhotonCandidate->SetLabelNegative(fESDArrayNeg[iGamma]);
2806                 continue;
2807                 }
2808                 fMCStackPos[iGamma] =  PhotonCandidate->GetMCLabelPositive();
2809                 fMCStackNeg[iGamma] =  PhotonCandidate->GetMCLabelNegative();
2810                 fESDArrayPos[iGamma] = PhotonCandidate->GetTrackLabelPositive();
2811                 fESDArrayNeg[iGamma] = PhotonCandidate->GetTrackLabelNegative();
2812
2813                 Bool_t AODLabelPos = kFALSE;
2814                 Bool_t AODLabelNeg = kFALSE;
2815
2816                 for(Int_t i = 0; i<fInputEvent->GetNumberOfTracks();i++){
2817                 AliAODTrack *tempDaughter = static_cast<AliAODTrack*>(fInputEvent->GetTrack(i));
2818                 if(!AODLabelPos){
2819                         if( tempDaughter->GetID() == PhotonCandidate->GetTrackLabelPositive() ){
2820                         PhotonCandidate->SetMCLabelPositive(abs(tempDaughter->GetLabel()));
2821                         PhotonCandidate->SetLabelPositive(i);
2822                         AODLabelPos = kTRUE;
2823                         }
2824                 }
2825                 if(!AODLabelNeg){
2826                         if( tempDaughter->GetID() == PhotonCandidate->GetTrackLabelNegative()){
2827                         PhotonCandidate->SetMCLabelNegative(abs(tempDaughter->GetLabel()));
2828                         PhotonCandidate->SetLabelNegative(i);
2829                         AODLabelNeg = kTRUE;
2830                         }
2831                 }
2832                 if(AODLabelNeg && AODLabelPos){
2833                         break;
2834                 }
2835                 }
2836                 if(!AODLabelPos || !AODLabelNeg){
2837                 cout<<"WARNING!!! AOD TRACKS NOT FOUND FOR"<<endl;
2838                 }
2839         }
2840
2841
2842         if(!mode){
2843                 delete[] fMCStackPos;
2844                 delete[] fMCStackNeg;
2845                 delete[] fESDArrayPos;
2846                 delete[] fESDArrayNeg;
2847         }
2848 }
2849
2850 void AliAnalysisTaskGammaConvCalo::SetLogBinningXTH2(TH2* histoRebin){
2851         TAxis *axisafter = histoRebin->GetXaxis();
2852         Int_t bins = axisafter->GetNbins();
2853         Double_t from = axisafter->GetXmin();
2854         Double_t to = axisafter->GetXmax();
2855         Double_t *newbins = new Double_t[bins+1];
2856         newbins[0] = from;
2857         Double_t factor = TMath::Power(to/from, 1./bins);
2858         for(Int_t i=1; i<=bins; ++i) newbins[i] = factor * newbins[i-1];
2859         axisafter->Set(bins, newbins);
2860         delete [] newbins;
2861 }
2862
2863 //________________________________________________________________________
2864 void AliAnalysisTaskGammaConvCalo::Terminate(const Option_t *)
2865 {
2866
2867   //fOutputContainer->Print(); // Will crash on GRID
2868 }
2869
2870 //________________________________________________________________________
2871 Int_t AliAnalysisTaskGammaConvCalo::GetSourceClassification(Int_t daughter, Int_t pdgCode){
2872
2873         if (daughter == 111) {
2874                 if (abs(pdgCode) == 310) return 1; // k0s
2875                 else if (abs(pdgCode) == 3122) return 2; // Lambda
2876                 else if (abs(pdgCode) == 130) return 3; // K0L
2877                 else if (abs(pdgCode) == 2212) return 4; // proton
2878                 else if (abs(pdgCode) == 2112) return 5; // neutron
2879                 else if (abs(pdgCode) == 211) return 6; // pion
2880                 else if (abs(pdgCode) == 321) return 7; // kaon
2881                 else if (abs(pdgCode) == 113 || abs(pdgCode) == 213 ) return 8; // rho 0,+,-
2882                 else if (abs(pdgCode) == 3222 || abs(pdgCode) == 3212 || abs(pdgCode) == 3112  ) return 9; // Sigma
2883                 else if (abs(pdgCode) == 2224 || abs(pdgCode) == 2214 || abs(pdgCode) == 2114 || abs(pdgCode) == 1114  ) return 10; // Delta
2884                 else if (abs(pdgCode) == 313 || abs(pdgCode) == 323   ) return 11; // K*
2885                 else return 15;
2886         }
2887         return 15;
2888
2889 }
2890
2891 // //________________________________________________________________________
2892 // Double_t AliAnalysisTaskGammaConvCalo::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
2893 // {
2894 //   // Get maximum energy of attached cell.
2895 //
2896 //   id = -1;
2897 //   Double_t maxe = 0;
2898 //   Int_t ncells = cluster->GetNCells();
2899 //   if (fEsdCells) {
2900 //     for (Int_t i=0; i<ncells; i++) {
2901 //       Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2902 //       if (e>maxe) {
2903 //         maxe = e;
2904 //         id   = cluster->GetCellAbsId(i);
2905 //       }
2906 //     }
2907 //   } else {
2908 //     for (Int_t i=0; i<ncells; i++) {
2909 //       Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2910 //       if (e>maxe)
2911 //         maxe = e;
2912 //       id   = cluster->GetCellAbsId(i);
2913 //     }
2914 //   }
2915 //   return maxe;
2916 // }