PWGGA/GammaConv/AliAnalysisTaskGammaConvV1.cxx

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