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