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