PWGGA/GammaConv/AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: Friederike Bock                                                                                    *
   5  * Version 1                                                              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //
  17
  18 #include <vector>
  19 #include "TParticle.h"
  20 #include "TPDGCode.h"
  21 #include "TMCProcess.h"
  22 #include "TDatabasePDG.h"
  23 #include "TList.h"
  24 #include "TChain.h"
  25 #include "TDirectory.h"
  26 #include "TTree.h"
  27 #include "TH1.h"
  28 #include "TH1F.h"
  29 #include "THnSparse.h"
  30 #include "TH2F.h"
  31 #include "AliStack.h"
  32 #include "AliAnalysisManager.h"
  33 #include "AliESDInputHandler.h"
  34 #include "AliESDtrack.h"
  35 #include "AliMCEvent.h"
  36 #include "AliStack.h"
  37 #include "AliMCEventHandler.h"
  38 #include "AliPID.h"
  39 #include "AliLog.h"
  40 #include "AliESDtrackCuts.h"
  41 #include "AliESDpidCuts.h"
  42 #include "AliMCEvent.h"
  43 #include "AliESDv0.h"
  44 #include "AliESDEvent.h"
  45 #include "AliESDpid.h"
  46 #include "AliKFParticle.h"
  47 #include "AliMCEventHandler.h"
  48 #include "AliKFVertex.h"
  49 #include "AliTriggerAnalysis.h"
  50 #include "AliCentrality.h"
  51 #include "AliMultiplicity.h"
  52 #include "AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero.h"
  53
  54
  55 ClassImp( AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero )
  56
  57 //-----------------------------------------------------------------------------------------------
  58 AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero():
  59         fV0Reader(NULL),
  60         fPionSelector(NULL),
  61         fBGHandler(NULL),
  62         fESDEvent(NULL),
  63         fMCEvent(NULL),
  64         fMCStack(NULL),
  65         fCutFolder(NULL),
  66         fESDList(NULL),
  67         fBackList(NULL),
  68         fMotherList(NULL),
  69         fTrueList(NULL),
  70         fMCList(NULL),
  71         fOutputContainer(0),
  72         fReaderGammas(NULL),
  73         fSelectorNegPionIndex(0),
  74         fSelectorPosPionIndex(0),
  75         fGoodConvGammas(NULL),
  76         fClusterCandidates(NULL),
  77         fNeutralPionCandidates(NULL),
  78         fGoodVirtualParticles(NULL),
  79         fEventCutArray(NULL),
  80         fGammaCutArray(NULL),
  81         fClusterCutArray(NULL),
  82         fPionCutArray(NULL),
  83         fNeutralPionMesonCutArray(NULL),
  84         fMesonCutArray(NULL),
  85         fEventCuts(NULL),
  86         fConversionCuts(NULL),
  87         fClusterCuts(NULL),
  88         fHistoConvGammaPt(NULL),
  89         fHistoConvGammaEta(NULL),
  90         fHistoClusterGammaPt(NULL),
  91         fHistoClusterGammaEta(NULL),
  92         fHistoNegPionPt(NULL),
  93         fHistoPosPionPt(NULL),
  94         fHistoNegPionPhi(NULL),
  95         fHistoPosPionPhi(NULL),
  96         fHistoNegPionEta(NULL),
  97         fHistoPosPionEta(NULL),
  98         fHistoNegPionClsTPC(NULL),
  99         fHistoPosPionClsTPC(NULL),
 100         fHistoPionDCAxy(NULL),
 101         fHistoPionDCAz(NULL),
 102         fHistoPionTPCdEdxNSigma(NULL),
 103         fHistoPionTPCdEdx(NULL),
 104         fHistoPionPionInvMassPt(NULL),
 105         fHistoGammaGammaInvMassPt(NULL),
 106         fHistoMotherInvMassPt(NULL),
 107         fTHnSparseMotherInvMassPtZM(NULL),
 108         fHistoMotherBackInvMassPt(NULL),
 109         fTHnSparseMotherBackInvMassPtZM(NULL),
 110         fHistoMCAllGammaPt(NULL),
 111         fHistoMCConvGammaPt(NULL),
 112         fHistoMCAllPosPionsPt(NULL),
 113         fHistoMCAllNegPionsPt(NULL),
 114         fHistoMCGammaFromNeutralMesonPt(NULL),
 115         fHistoMCPosPionsFromNeutralMesonPt(NULL),
 116         fHistoMCNegPionsFromNeutralMesonPt(NULL),
 117         fHistoMCEtaPiPlPiMiPiZeroPt(NULL),
 118         fHistoMCEtaPiPlPiMiPiZeroInAccPt(NULL),
 119         fHistoMCOmegaPiPlPiMiPiZeroPt(NULL),
 120         fHistoMCOmegaPiPlPiMiPiZeroInAccPt(NULL),
 121         fHistoTrueMotherPiPlPiMiPiZeroInvMassPt(NULL),
 122         fHistoTrueMotherGammaGammaInvMassPt(NULL),
 123         fHistoTrueMotherGammaGammaFromEtaInvMassPt(NULL),
 124         fHistoTrueMotherGammaGammaFromOmegaInvMassPt(NULL),
 125         fHistoTrueConvGammaPt(NULL),
 126         fHistoTrueConvGammaFromNeutralMesonPt(NULL),
 127         fHistoTrueClusterGammaPt(NULL),
 128         fHistoTrueClusterGammaFromNeutralMesonPt(NULL),
 129         fHistoTruePosPionPt(NULL),
 130         fHistoTruePosPionFromNeutralMesonPt(NULL),
 131         fHistoTrueNegPionPt(NULL),
 132         fHistoTrueNegPionFromNeutralMesonPt(NULL),
 133         fHistoTruePionPionInvMassPt(NULL),
 134         fHistoTruePionPionFromSameMotherInvMassPt(NULL),
 135         fHistoTruePionPionFromEtaInvMassPt(NULL),
 136         fHistoTruePionPionFromOmegaInvMassPt(NULL),
 137         fHistoNEvents(NULL),
 138         fHistoNGoodESDTracks(NULL),
 139         fProfileEtaShift(NULL),
 140         fRandom(0),
 141         fnCuts(0),
 142         fiCut(0),
 143         fNumberOfESDTracks(0),
 144         fMoveParticleAccordingToVertex(kFALSE),
 145         fIsHeavyIon(0),
 146         fDoMesonAnalysis(kTRUE),
 147         fDoMesonQA(kFALSE),
 148         fIsFromMBHeader(kTRUE),
 149         fIsMC(kFALSE),
 150         fNeutralPionMode(0)
 151 {
 152
 153 }
 154
 155 //-----------------------------------------------------------------------------------------------
 156 AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero( const char* name ):
 157         AliAnalysisTaskSE(name),
 158         fV0Reader(NULL),
 159         fPionSelector(NULL),
 160         fBGHandler(NULL),
 161         fESDEvent(NULL),
 162         fMCEvent(NULL),
 163         fMCStack(NULL),
 164         fCutFolder(NULL),
 165         fESDList(NULL),
 166         fBackList(NULL),
 167         fMotherList(NULL),
 168         fTrueList(NULL),
 169         fMCList(NULL),
 170         fOutputContainer(0),
 171         fReaderGammas(NULL),
 172         fSelectorNegPionIndex(0),
 173         fSelectorPosPionIndex(0),
 174         fGoodConvGammas(NULL),
 175         fClusterCandidates(NULL),
 176         fNeutralPionCandidates(NULL),
 177         fGoodVirtualParticles(NULL),
 178         fEventCutArray(NULL),
 179         fGammaCutArray(NULL),
 180         fClusterCutArray(NULL),
 181         fPionCutArray(NULL),
 182         fNeutralPionMesonCutArray(NULL),
 183         fMesonCutArray(NULL),
 184         fEventCuts(NULL),
 185         fConversionCuts(NULL),
 186         fClusterCuts(NULL),
 187         fHistoConvGammaPt(NULL),
 188         fHistoConvGammaEta(NULL),
 189         fHistoClusterGammaPt(NULL),
 190         fHistoClusterGammaEta(NULL),
 191         fHistoNegPionPt(NULL),
 192         fHistoPosPionPt(NULL),
 193         fHistoNegPionPhi(NULL),
 194         fHistoPosPionPhi(NULL),
 195         fHistoNegPionEta(NULL),
 196         fHistoPosPionEta(NULL),
 197         fHistoNegPionClsTPC(NULL),
 198         fHistoPosPionClsTPC(NULL),
 199         fHistoPionDCAxy(NULL),
 200         fHistoPionDCAz(NULL),
 201         fHistoPionTPCdEdxNSigma(NULL),
 202         fHistoPionTPCdEdx(NULL),
 203         fHistoPionPionInvMassPt(NULL),
 204         fHistoGammaGammaInvMassPt(NULL),
 205         fHistoMotherInvMassPt(NULL),
 206         fTHnSparseMotherInvMassPtZM(NULL),
 207         fHistoMotherBackInvMassPt(NULL),
 208         fTHnSparseMotherBackInvMassPtZM(NULL),
 209         fHistoMCAllGammaPt(NULL),
 210         fHistoMCConvGammaPt(NULL),
 211         fHistoMCAllPosPionsPt(NULL),
 212         fHistoMCAllNegPionsPt(NULL),
 213         fHistoMCGammaFromNeutralMesonPt(NULL),
 214         fHistoMCPosPionsFromNeutralMesonPt(NULL),
 215         fHistoMCNegPionsFromNeutralMesonPt(NULL),
 216         fHistoMCEtaPiPlPiMiPiZeroPt(NULL),
 217         fHistoMCEtaPiPlPiMiPiZeroInAccPt(NULL),
 218         fHistoMCOmegaPiPlPiMiPiZeroPt(NULL),
 219         fHistoMCOmegaPiPlPiMiPiZeroInAccPt(NULL),
 220         fHistoTrueMotherPiPlPiMiPiZeroInvMassPt(NULL),
 221         fHistoTrueMotherGammaGammaInvMassPt(NULL),
 222         fHistoTrueMotherGammaGammaFromEtaInvMassPt(NULL),
 223         fHistoTrueMotherGammaGammaFromOmegaInvMassPt(NULL),
 224         fHistoTrueConvGammaPt(NULL),
 225         fHistoTrueConvGammaFromNeutralMesonPt(NULL),
 226         fHistoTrueClusterGammaPt(NULL),
 227         fHistoTrueClusterGammaFromNeutralMesonPt(NULL),
 228         fHistoTruePosPionPt(NULL),
 229         fHistoTruePosPionFromNeutralMesonPt(NULL),
 230         fHistoTrueNegPionPt(NULL),
 231         fHistoTrueNegPionFromNeutralMesonPt(NULL),
 232         fHistoTruePionPionInvMassPt(NULL),
 233         fHistoTruePionPionFromSameMotherInvMassPt(NULL),
 234         fHistoTruePionPionFromEtaInvMassPt(NULL),
 235         fHistoTruePionPionFromOmegaInvMassPt(NULL),
 236         fHistoNEvents(NULL),
 237         fHistoNGoodESDTracks(NULL),
 238         fProfileEtaShift(NULL),
 239         fRandom(0),
 240         fnCuts(0),
 241         fiCut(0),
 242         fNumberOfESDTracks(0),
 243         fMoveParticleAccordingToVertex(kFALSE),
 244         fIsHeavyIon(0),
 245         fDoMesonAnalysis(kTRUE),
 246         fDoMesonQA(kFALSE),
 247         fIsFromMBHeader(kTRUE),
 248         fIsMC(kFALSE),
 249         fNeutralPionMode(0)
 250 {
 251    DefineOutput(1, TList::Class());
 252 }
 253
 254 //-----------------------------------------------------------------------------------------------
 255 AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::~AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero()
 256 {
 257         //
 258         // virtual destructor
 259         //
 260         cout<<"Destructor"<<endl;
 261         if(fGoodConvGammas){
 262                 delete fGoodConvGammas;
 263                 fGoodConvGammas = 0x0;
 264         }
 265         if(fClusterCandidates){
 266                 delete fClusterCandidates;
 267                 fClusterCandidates = 0x0;
 268         }
 269
 270         if(fNeutralPionCandidates){
 271                 delete fNeutralPionCandidates;
 272                 fNeutralPionCandidates = 0x0;
 273         }
 274
 275         if(fGoodVirtualParticles){
 276                 delete fGoodVirtualParticles;
 277                 fGoodVirtualParticles = 0x0;
 278         }
 279         if(fBGHandler){
 280                 delete[] fBGHandler;
 281                 fBGHandler = 0x0;
 282         }
 283 }
 284 //___________________________________________________________
 285 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::InitBack(){
 286
 287         const Int_t nDim = 4;
 288         Int_t nBins[nDim] = {500,250,7,4};
 289         Double_t xMin[nDim] = {0.4,0, 0,0};
 290         Double_t xMax[nDim] = {0.9,25,7,4};
 291
 292         fTHnSparseMotherInvMassPtZM = new THnSparseF*[fnCuts];
 293         fTHnSparseMotherBackInvMassPtZM = new THnSparseF*[fnCuts];
 294
 295         fBGHandler = new AliGammaConversionAODBGHandler*[fnCuts];
 296
 297         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 298
 299                 TString cutstringEvent          = ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber();
 300                 TString cutstringPion           = ((AliPrimaryPionCuts*)fPionCutArray->At(iCut))->GetCutNumber();
 301                 TString cutstringConvGamma = "";
 302                 if (fNeutralPionMode < 2)  cutstringConvGamma = ((AliConversionPhotonCuts*)fGammaCutArray->At(iCut))->GetCutNumber();
 303                 TString cutstringCaloGamma = "";
 304                 if (fNeutralPionMode > 0)  cutstringCaloGamma = ((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutNumber();
 305                 TString cutstringNeutralPion= ((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(iCut))->GetCutNumber();
 306                 TString cutstringMeson          = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
 307
 308                 TString fullCutString = "";
 309                 if (fNeutralPionMode == 0) fullCutString = Form("%i_%s_%s_%s_%s_%s",fNeutralPionMode,cutstringEvent.Data(),cutstringConvGamma.Data(),cutstringNeutralPion.Data(), cutstringPion.Data(),cutstringMeson.Data());
 310                 else if (fNeutralPionMode == 1) fullCutString = Form("%i_%s_%s_%s_%s_%s_%s",fNeutralPionMode,cutstringEvent.Data(),cutstringConvGamma.Data(),cutstringCaloGamma.Data(),cutstringNeutralPion.Data(), cutstringPion.Data(),cutstringMeson.Data());
 311                 else if (fNeutralPionMode == 2) fullCutString = Form("%i_%s_%s_%s_%s_%s",fNeutralPionMode,cutstringEvent.Data(),cutstringCaloGamma.Data(),cutstringNeutralPion.Data(), cutstringPion.Data(),cutstringMeson.Data());
 312
 313                 TString nameBackList = Form("%s Back histograms",fullCutString.Data());
 314                 TString nameMotherList = Form("%s Mother histograms",fullCutString.Data());
 315
 316                 Int_t collisionSystem = atoi((TString)(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber())(0,1));
 317                 Int_t centMin = atoi((TString)(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber())(1,1));
 318                 Int_t centMax = atoi((TString)(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber())(2,1));
 319
 320                 if(collisionSystem == 1 || collisionSystem == 2 ||
 321                         collisionSystem == 5 || collisionSystem == 8 ||
 322                         collisionSystem == 9){
 323                         centMin = centMin*10;
 324                         centMax = centMax*10;
 325                 }
 326                 else if(collisionSystem == 3 || collisionSystem == 6){
 327                         centMin = centMin*5;
 328                         centMax = centMax*5;
 329                 }
 330                 else if(collisionSystem == 4 || collisionSystem == 7){
 331                         centMin = ((centMin*5)+45);
 332                         centMax = ((centMax*5)+45);
 333                 }
 334
 335
 336                 fBackList[iCut] = new TList();
 337
 338                 fBackList[iCut]->SetName(nameBackList.Data());
 339                 fBackList[iCut]->SetOwner(kTRUE);
 340                 fCutFolder[iCut]->Add(fBackList[iCut]);
 341
 342                 fTHnSparseMotherBackInvMassPtZM[iCut] = new THnSparseF("Back_Back_InvMass_Pt_z_m","Back_Back_InvMass_Pt_z_m",nDim,nBins,xMin,xMax);
 343                 fBackList[iCut]->Add(fTHnSparseMotherBackInvMassPtZM[iCut]);
 344
 345                 fMotherList[iCut] = new TList();
 346                 fMotherList[iCut]->SetName(nameMotherList.Data());
 347                 fMotherList[iCut]->SetOwner(kTRUE);
 348                 fCutFolder[iCut]->Add(fMotherList[iCut]);
 349
 350                 fTHnSparseMotherInvMassPtZM[iCut] = new THnSparseF("Back_Mother_InvMass_Pt_z_m","Back_Mother_InvMass_Pt_z_m",nDim,nBins,xMin,xMax);
 351                 fMotherList[iCut]->Add(fTHnSparseMotherInvMassPtZM[iCut]);
 352
 353
 354                 fBGHandler[iCut] = new AliGammaConversionAODBGHandler(  collisionSystem,centMin,centMax,
 355                                                                                                                                 ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->GetNumberOfBGEvents(),
 356                                                                                                                                 ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity());
 357         }
 358 }
 359
 360 //______________________________________________________________________
 361 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::UserCreateOutputObjects()
 362 {
 363         //
 364         // Create ouput objects
 365         //
 366
 367         // Create the output container
 368         if(fOutputContainer != NULL){
 369                 delete fOutputContainer;
 370                 fOutputContainer = NULL;
 371         }
 372         if(fOutputContainer == NULL){
 373                 fOutputContainer = new TList();
 374                 fOutputContainer->SetOwner(kTRUE);
 375         }
 376
 377         fGoodConvGammas = new TList();
 378         fClusterCandidates = new TList();
 379         fNeutralPionCandidates = new TList();
 380         fGoodVirtualParticles = new TList();
 381
 382         fCutFolder                              = new TList*[fnCuts];
 383         fESDList                                = new TList*[fnCuts];
 384         fBackList                               = new TList*[fnCuts];
 385         fMotherList                     = new TList*[fnCuts];
 386         fHistoNEvents                   = new TH1I*[fnCuts];
 387         fHistoNGoodESDTracks    = new TH1I*[fnCuts];
 388         fProfileEtaShift                = new TProfile*[fnCuts];
 389         if (fNeutralPionMode < 2){
 390                 fHistoConvGammaPt               = new TH1F*[fnCuts];
 391                 fHistoConvGammaEta              = new TH1F*[fnCuts];
 392         }
 393         if (fNeutralPionMode > 0){
 394                 fHistoClusterGammaPt            = new TH1F*[fnCuts];
 395                 fHistoClusterGammaEta           = new TH1F*[fnCuts];
 396         }
 397         fHistoNegPionPt                 = new TH1F*[fnCuts];
 398         fHistoPosPionPt                 = new TH1F*[fnCuts];
 399         fHistoNegPionPhi                = new TH1F*[fnCuts];
 400         fHistoPosPionPhi                = new TH1F*[fnCuts];
 401         fHistoPionPionInvMassPt = new TH2F*[fnCuts];
 402
 403         if( fDoMesonQA ) {
 404                 fHistoNegPionEta                = new TH1F*[fnCuts];
 405                 fHistoPosPionEta                = new TH1F*[fnCuts];
 406                 fHistoNegPionClsTPC             = new TH2F*[fnCuts];
 407                 fHistoPosPionClsTPC             = new TH2F*[fnCuts];
 408                 fHistoPionDCAxy                 = new TH2F*[fnCuts];
 409                 fHistoPionDCAz                  = new TH2F*[fnCuts];
 410                 fHistoPionTPCdEdxNSigma = new TH2F*[fnCuts];
 411                 fHistoPionTPCdEdx               = new TH2F*[fnCuts];
 412         }
 413         fHistoGammaGammaInvMassPt       = new TH2F*[fnCuts];
 414         fHistoMotherInvMassPt           = new TH2F*[fnCuts];
 415         fHistoMotherBackInvMassPt       = new TH2F*[fnCuts];
 416
 417         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 418                 TString cutstringEvent          = ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber();
 419                 TString cutstringPion           = ((AliPrimaryPionCuts*)fPionCutArray->At(iCut))->GetCutNumber();
 420                 TString cutstringConvGamma = "";
 421                 if (fNeutralPionMode < 2)  cutstringConvGamma = ((AliConversionPhotonCuts*)fGammaCutArray->At(iCut))->GetCutNumber();
 422                 TString cutstringCaloGamma = "";
 423                 if (fNeutralPionMode > 0)  cutstringCaloGamma = ((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutNumber();
 424                 TString cutstringNeutralPion= ((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(iCut))->GetCutNumber();
 425                 TString cutstringMeson          = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
 426
 427                 TString fullCutString = "";
 428                 if (fNeutralPionMode == 0) fullCutString = Form("%i_%s_%s_%s_%s_%s",fNeutralPionMode,cutstringEvent.Data(),cutstringConvGamma.Data(),cutstringNeutralPion.Data(), cutstringPion.Data(),cutstringMeson.Data());
 429                 else if (fNeutralPionMode == 1) fullCutString = Form("%i_%s_%s_%s_%s_%s_%s",fNeutralPionMode,cutstringEvent.Data(),cutstringConvGamma.Data(),cutstringCaloGamma.Data(),cutstringNeutralPion.Data(), cutstringPion.Data(),cutstringMeson.Data());
 430                 else if (fNeutralPionMode == 2) fullCutString = Form("%i_%s_%s_%s_%s_%s",fNeutralPionMode,cutstringEvent.Data(),cutstringCaloGamma.Data(),cutstringNeutralPion.Data(), cutstringPion.Data(),cutstringMeson.Data());
 431                 TString nameCutFolder = Form("Cut Number %s", fullCutString.Data());
 432                 TString nameESDList = Form("%s ESD histograms", fullCutString.Data());
 433
 434
 435                 fCutFolder[iCut] = new TList();
 436                 fCutFolder[iCut]->SetName(nameCutFolder.Data());
 437                 fCutFolder[iCut]->SetOwner(kTRUE);
 438                 fOutputContainer->Add(fCutFolder[iCut]);
 439
 440                 fESDList[iCut] = new TList();
 441                 fESDList[iCut]->SetName(nameESDList.Data());
 442                 fESDList[iCut]->SetOwner(kTRUE);
 443
 444                 fHistoNEvents[iCut] = new TH1I("NEvents","NEvents",9,-0.5,8.5);
 445                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(1,"Accepted");
 446                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(2,"Centrality");
 447                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(3,"Missing MC");
 448                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(4,"Trigger");
 449                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(5,"Vertex Z");
 450                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(6,"Cont. Vertex");
 451                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(7,"Pile-Up");
 452                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(8,"no SDD");
 453                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(9,"no V0AND");
 454                 fESDList[iCut]->Add(fHistoNEvents[iCut]);
 455
 456                 if(fIsHeavyIon>0) fHistoNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",3000,0,3000);
 457                         else fHistoNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",200,0,200);
 458                 fESDList[iCut]->Add(fHistoNGoodESDTracks[iCut]);
 459
 460                 fProfileEtaShift[iCut] = new TProfile("Eta Shift","Eta Shift",1, -0.5,0.5);
 461                 fESDList[iCut]->Add(fProfileEtaShift[iCut]);
 462                 if (fNeutralPionMode < 2){
 463                         fHistoConvGammaPt[iCut] = new TH1F("ESD_ConvGamma_Pt","ESD_ConvGamma_Pt",250,0,25);
 464                         fESDList[iCut]->Add(fHistoConvGammaPt[iCut]);
 465                         fHistoConvGammaEta[iCut] = new TH1F("ESD_ConvGamma_Eta","ESD_ConvGamma_Eta",600,-1.5,1.5);
 466                         fESDList[iCut]->Add(fHistoConvGammaEta[iCut]);
 467                 }
 468                 if (fNeutralPionMode > 0){
 469                         fHistoClusterGammaPt[iCut] = new TH1F("ESD_ClusterGamma_Pt","ESD_ClusterGamma_Pt",250,0,25);
 470                         fESDList[iCut]->Add(fHistoClusterGammaPt[iCut]);
 471                         fHistoClusterGammaEta[iCut] = new TH1F("ESD_ClusterGamma_Eta","ESD_ClusterGamma_Eta",600,-1.5,1.5);
 472                         fESDList[iCut]->Add(fHistoClusterGammaEta[iCut]);
 473                 }
 474                 fHistoNegPionPt[iCut] = new TH1F("ESD_PrimaryNegPions_Pt","ESD_PrimaryNegPions_Pt",1000,0,25);
 475                 fESDList[iCut]->Add(fHistoNegPionPt[iCut]);
 476                 fHistoPosPionPt[iCut] = new TH1F("ESD_PrimaryPosPions_Pt","ESD_PrimaryPosPions_Pt",1000,0,25);
 477                 fESDList[iCut]->Add(fHistoPosPionPt[iCut]);
 478                 fHistoNegPionPhi[iCut] = new TH1F("ESD_PrimaryNegPions_Phi","ESD_PrimaryNegPions_Phi",360,0,2*TMath::Pi());
 479                 fESDList[iCut]->Add(fHistoNegPionPhi[iCut]);
 480                 fHistoPosPionPhi[iCut] = new TH1F("ESD_PrimaryPosPions_Phi","ESD_PrimaryPosPions_Phi",360,0,2*TMath::Pi());
 481                 fESDList[iCut]->Add(fHistoPosPionPhi[iCut]);
 482                 fHistoPionPionInvMassPt[iCut] = new TH2F("ESD_PiPlusPiNeg_InvMassPt","ESD_PiPlusPiNeg_InvMassPt",2000,0.,2.,200,0.,20.);
 483                 fESDList[iCut]->Add(fHistoPionPionInvMassPt[iCut]);
 484
 485                 if ( fDoMesonQA ) {
 486                         fHistoNegPionEta[iCut] = new TH1F("ESD_PrimaryNegPions_Eta","ESD_PrimaryNegPions_Eta",600,-1.5,1.5);
 487                         fESDList[iCut]->Add(fHistoNegPionEta[iCut]);
 488                         fHistoPosPionEta[iCut] = new TH1F("ESD_PrimaryPosPions_Eta","ESD_PrimaryPosPions_Eta",600,-1.5,1.5);
 489                         fESDList[iCut]->Add(fHistoPosPionEta[iCut]);
 490                         fHistoNegPionClsTPC[iCut]  = new TH2F("ESD_PrimaryNegPions_ClsTPC","ESD_PrimaryNegPions_ClsTPC",100,0,1,400,0.,10.);
 491                         fESDList[iCut]->Add(fHistoNegPionClsTPC[iCut]);
 492                         fHistoPosPionClsTPC[iCut]  = new TH2F("ESD_PrimaryPosPions_ClsTPC","ESD_PrimaryPosPions_ClsTPC",100,0,1,400,0.,10.);
 493                         fESDList[iCut]->Add(fHistoPosPionClsTPC[iCut]);
 494                         fHistoPionDCAxy[iCut] = new TH2F("ESD_PrimaryPions_DCAxy","ESD_PrimaryPions_DCAxy",800,-4.0,4.0,400,0.,10.);
 495                         fESDList[iCut]->Add(fHistoPionDCAxy[iCut]);
 496                         fHistoPionDCAz[iCut]  = new TH2F("ESD_PrimaryPions_DCAz","ESD_PrimaryPions_DCAz",800,-4.0,4.0,400,0.,10.);
 497                         fESDList[iCut]->Add(fHistoPionDCAz[iCut]);
 498                         fHistoPionTPCdEdxNSigma[iCut] = new TH2F("ESD_PrimaryPions_TPCdEdx","ESD_PrimaryPions_TPCdEdx",150,0.05,20,400,-10,10);
 499                         fESDList[iCut]->Add(fHistoPionTPCdEdxNSigma[iCut]);
 500                         fHistoPionTPCdEdx[iCut] =new TH2F("ESD_PrimaryPions_TPCdEdxSignal","ESD_PrimaryPions_TPCdEdxSignal" ,150,0.05,20.0,800,0.0,200);
 501                         fESDList[iCut]->Add(fHistoPionTPCdEdx[iCut]);
 502                 }
 503                 fHistoGammaGammaInvMassPt[iCut] = new TH2F("ESD_GammaGamma_InvMass_Pt","ESD_GammaGamma_InvMass_Pt",450,0.,0.45,250,0,25);
 504                 fESDList[iCut]->Add(fHistoGammaGammaInvMassPt[iCut]);
 505                 fHistoMotherInvMassPt[iCut] = new TH2F("ESD_Mother_InvMass_Pt","ESD_Mother_InvMass_Pt",500,0.4,0.9,250,0,25);
 506                 fESDList[iCut]->Add(fHistoMotherInvMassPt[iCut]);
 507                 fHistoMotherBackInvMassPt[iCut] = new TH2F("ESD_Background_InvMass_Pt","ESD_Background_InvMass_Pt",500,0.4,0.9,250,0,25);
 508                 fESDList[iCut]->Add(fHistoMotherBackInvMassPt[iCut]);
 509
 510                 if ( fDoMesonQA ) {
 511                         TAxis *AxisAfter = fHistoPionTPCdEdxNSigma[iCut]->GetXaxis();
 512                         Int_t bins = AxisAfter->GetNbins();
 513                         Double_t from = AxisAfter->GetXmin();
 514                         Double_t to = AxisAfter->GetXmax();
 515                         Double_t *newBins = new Double_t[bins+1];
 516                         newBins[0] = from;
 517                         Double_t factor = TMath::Power(to/from, 1./bins);
 518                         for(Int_t i=1; i<=bins; ++i) newBins[i] = factor * newBins[i-1];
 519
 520                         AxisAfter->Set(bins, newBins);
 521                         AxisAfter = fHistoPionTPCdEdx[iCut]->GetXaxis();
 522                         AxisAfter->Set(bins, newBins);
 523
 524                         delete [] newBins;
 525                 }
 526
 527                 fCutFolder[iCut]->Add(fESDList[iCut]);
 528
 529         }
 530
 531         if( fIsMC ){
 532                 // MC Histogramms
 533                 fMCList = new TList*[fnCuts];
 534                 // True Histogramms
 535                 fTrueList = new TList*[fnCuts];
 536                 if (fNeutralPionMode < 2){
 537                         fHistoTrueConvGammaPt = new TH1F*[fnCuts];
 538                         fHistoTrueConvGammaFromNeutralMesonPt = new TH1F*[fnCuts];
 539                 }
 540                 if (fNeutralPionMode > 0){
 541                         fHistoTrueClusterGammaPt = new TH1F*[fnCuts];
 542                         fHistoTrueClusterGammaFromNeutralMesonPt = new TH1F*[fnCuts];
 543                 }
 544                 fHistoTruePosPionPt  = new TH1F*[fnCuts];
 545                 fHistoTrueNegPionPt  = new TH1F*[fnCuts];
 546                 fHistoTruePosPionFromNeutralMesonPt  = new TH1F*[fnCuts];
 547                 fHistoTrueNegPionFromNeutralMesonPt  = new TH1F*[fnCuts];
 548
 549
 550                 fHistoMCAllGammaPt  = new TH1F*[fnCuts];
 551                 if (fNeutralPionMode < 2){
 552                         fHistoMCConvGammaPt = new TH1F*[fnCuts];
 553                 }
 554                 fHistoMCAllPosPionsPt = new TH1F*[fnCuts];
 555                 fHistoMCAllNegPionsPt = new TH1F*[fnCuts];
 556                 fHistoMCGammaFromNeutralMesonPt  = new TH1F*[fnCuts];
 557                 fHistoMCPosPionsFromNeutralMesonPt = new TH1F*[fnCuts];
 558                 fHistoMCNegPionsFromNeutralMesonPt = new TH1F*[fnCuts];
 559
 560 //              hMCPi0DalitzGammaPt    = new TH1F*[fnCuts];
 561 //              hMCPi0DalitzElectronPt = new TH1F*[fnCuts];
 562 //              hMCPi0DalitzPositronPt = new TH1F*[fnCuts];
 563
 564                 fHistoMCEtaPiPlPiMiPiZeroPt = new TH1F*[fnCuts];
 565                 fHistoMCEtaPiPlPiMiPiZeroInAccPt = new TH1F*[fnCuts];
 566                 fHistoMCOmegaPiPlPiMiPiZeroPt = new TH1F*[fnCuts];
 567                 fHistoMCOmegaPiPlPiMiPiZeroInAccPt = new TH1F*[fnCuts];
 568
 569                 fHistoTrueMotherPiPlPiMiPiZeroInvMassPt = new TH2F*[fnCuts];
 570                 fHistoTrueMotherGammaGammaInvMassPt = new TH2F*[fnCuts];
 571                 fHistoTrueMotherGammaGammaFromEtaInvMassPt = new TH2F*[fnCuts];
 572                 fHistoTrueMotherGammaGammaFromOmegaInvMassPt = new TH2F*[fnCuts];
 573 //              if (fDoMesonQA){
 574                         fHistoTruePionPionInvMassPt =                   new TH2F*[fnCuts];
 575                         fHistoTruePionPionFromSameMotherInvMassPt =     new TH2F*[fnCuts];
 576                         fHistoTruePionPionFromEtaInvMassPt =    new TH2F*[fnCuts];
 577                         fHistoTruePionPionFromOmegaInvMassPt =  new TH2F*[fnCuts];
 578 //              }
 579
 580                 for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 581                         TString cutstringEvent          = ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber();
 582                         TString cutstringPion           = ((AliPrimaryPionCuts*)fPionCutArray->At(iCut))->GetCutNumber();
 583                         TString cutstringConvGamma = "";
 584                         if (fNeutralPionMode < 2)  cutstringConvGamma = ((AliConversionPhotonCuts*)fGammaCutArray->At(iCut))->GetCutNumber();
 585                         TString cutstringCaloGamma = "";
 586                         if (fNeutralPionMode > 0)  cutstringCaloGamma = ((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutNumber();
 587                         TString cutstringNeutralPion= ((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(iCut))->GetCutNumber();
 588                         TString cutstringMeson          = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
 589
 590                         TString fullCutString = "";
 591                         if (fNeutralPionMode == 0) fullCutString = Form("%i_%s_%s_%s_%s_%s",fNeutralPionMode,cutstringEvent.Data(),cutstringConvGamma.Data(),cutstringNeutralPion.Data(), cutstringPion.Data(),cutstringMeson.Data());
 592                         else if (fNeutralPionMode == 1) fullCutString = Form("%i_%s_%s_%s_%s_%s_%s",fNeutralPionMode,cutstringEvent.Data(),cutstringConvGamma.Data(),cutstringCaloGamma.Data(),cutstringNeutralPion.Data(), cutstringPion.Data(),cutstringMeson.Data());
 593                         else if (fNeutralPionMode == 2) fullCutString = Form("%i_%s_%s_%s_%s_%s",fNeutralPionMode,cutstringEvent.Data(),cutstringCaloGamma.Data(),cutstringNeutralPion.Data(), cutstringPion.Data(),cutstringMeson.Data());
 594                         TString nameMCList = Form("%s MC histograms", fullCutString.Data());
 595                         TString nameTrueRecList = Form("%s True histograms", fullCutString.Data());
 596
 597                         fMCList[iCut] = new TList();
 598                         fMCList[iCut]->SetName(nameMCList.Data());
 599                         fMCList[iCut]->SetOwner(kTRUE);
 600                         fCutFolder[iCut]->Add(fMCList[iCut]);
 601
 602                         fHistoMCAllGammaPt[iCut] = new TH1F("MC_AllGamma_Pt","MC_AllGamma_Pt",250,0,25);
 603                         fMCList[iCut]->Add(fHistoMCAllGammaPt[iCut]);
 604                         if (fNeutralPionMode < 2){
 605                                 fHistoMCConvGammaPt[iCut] = new TH1F("MC_ConvGamma_Pt","MC_ConvGamma_Pt",250,0,25);
 606                                 fMCList[iCut]->Add(fHistoMCConvGammaPt[iCut]);
 607                         }
 608
 609                         fHistoMCAllPosPionsPt[iCut] = new TH1F("MC_AllPosPions_Pt","MC_AllPosPions_Pt",1000,0,25);
 610                         fMCList[iCut]->Add(fHistoMCAllPosPionsPt[iCut]);
 611                         fHistoMCAllNegPionsPt[iCut] = new TH1F("MC_AllNegPions_Pt","MC_AllNegPions_Pt",1000,0,25);
 612                         fMCList[iCut]->Add(fHistoMCAllNegPionsPt[iCut]);
 613                         fHistoMCGammaFromNeutralMesonPt[iCut] = new TH1F("MC_GammaFromNeutralMeson_Pt","MC_GammaFromNeutralMeson_Pt",250,0,25);
 614                         fMCList[iCut]->Add(fHistoMCGammaFromNeutralMesonPt[iCut]);
 615                         fHistoMCPosPionsFromNeutralMesonPt[iCut] = new TH1F("MC_PosPionsFromNeutralMeson_Pt","MC_PosPionsFromNeutralMeson_Pt",1000,0,25);
 616                         fMCList[iCut]->Add(fHistoMCPosPionsFromNeutralMesonPt[iCut]);
 617                         fHistoMCNegPionsFromNeutralMesonPt[iCut] = new TH1F("MC_NegPionsFromNeutralMeson_Pt","MC_NegPionsFromNeutralMeson_Pt",1000,0,25);
 618                         fMCList[iCut]->Add(fHistoMCNegPionsFromNeutralMesonPt[iCut]);
 619
 620                         fHistoMCEtaPiPlPiMiPiZeroPt[iCut] = new TH1F("MC_Eta_Pt","MC_Eta_Pt",250,0,25);
 621                         fHistoMCEtaPiPlPiMiPiZeroPt[iCut]->Sumw2();
 622                         fMCList[iCut]->Add(fHistoMCEtaPiPlPiMiPiZeroPt[iCut]);
 623
 624                         fHistoMCEtaPiPlPiMiPiZeroInAccPt[iCut] = new TH1F("MC_EtaInAcc_Pt","MC_EtaInAcc_Pt",250,0,25);
 625                         fHistoMCEtaPiPlPiMiPiZeroInAccPt[iCut]->Sumw2();
 626                         fMCList[iCut]->Add(fHistoMCEtaPiPlPiMiPiZeroInAccPt[iCut]);
 627
 628                         fHistoMCOmegaPiPlPiMiPiZeroPt[iCut] = new TH1F("MC_Omega_Pt","MC_Omega_Pt",250,0,25);
 629                         fHistoMCOmegaPiPlPiMiPiZeroPt[iCut]->Sumw2();
 630                         fMCList[iCut]->Add(fHistoMCOmegaPiPlPiMiPiZeroPt[iCut]);
 631
 632                         fHistoMCOmegaPiPlPiMiPiZeroInAccPt[iCut] = new TH1F("MC_OmegaInAcc_Pt","MC_OmegaInAcc_Pt",250,0,25);
 633                         fHistoMCOmegaPiPlPiMiPiZeroInAccPt[iCut]->Sumw2();
 634                         fMCList[iCut]->Add(fHistoMCOmegaPiPlPiMiPiZeroInAccPt[iCut]);
 635
 636                         fTrueList[iCut] = new TList();
 637                         fTrueList[iCut]->SetName(nameTrueRecList.Data());
 638                         fTrueList[iCut]->SetOwner(kTRUE);
 639                         fCutFolder[iCut]->Add(fTrueList[iCut]);
 640
 641                         if (fNeutralPionMode < 2){
 642                                 fHistoTrueConvGammaPt[iCut] = new TH1F("ESD_TrueConvGamma_Pt","ESD_TrueConvGamma_Pt",250,0,25);
 643                                 fTrueList[iCut]->Add(fHistoTrueConvGammaPt[iCut]);
 644                                 fHistoTrueConvGammaFromNeutralMesonPt[iCut] = new TH1F("ESD_TrueConvGammaFromNeutralMeson_Pt","ESD_TrueConvGammaFromNeutralMeson_Pt",250,0,25);
 645                                 fTrueList[iCut]->Add(fHistoTrueConvGammaFromNeutralMesonPt[iCut]);
 646                         }
 647                         if (fNeutralPionMode > 0){
 648                                 fHistoTrueClusterGammaPt[iCut] = new TH1F("ESD_TrueClusterGamma_Pt","ESD_TrueClusterGamma_Pt",250,0,25);
 649                                 fTrueList[iCut]->Add(fHistoTrueClusterGammaPt[iCut]);
 650                                 fHistoTrueClusterGammaFromNeutralMesonPt[iCut] = new TH1F("ESD_TrueClusterGammaFromNeutralMeson_Pt","ESD_TrueClusterGammaFromNeutralMeson_Pt",250,0,25);
 651                                 fTrueList[iCut]->Add(fHistoTrueClusterGammaFromNeutralMesonPt[iCut]);
 652                         }
 653                         fHistoTruePosPionPt[iCut] = new TH1F("ESD_TruePosPion_Pt","ESD_TruePosPion_Pt",1000,0,25);
 654                         fTrueList[iCut]->Add(fHistoTruePosPionPt[iCut]);
 655                         fHistoTrueNegPionPt[iCut] = new TH1F("ESD_TrueNegPion_Pt","ESD_TrueNegPion_Pt",1000,0,25);
 656                         fTrueList[iCut]->Add(fHistoTrueNegPionPt[iCut]);
 657
 658                         fHistoTrueNegPionFromNeutralMesonPt[iCut] = new TH1F("ESD_TrueNegPionFromNeutralMeson_Pt","ESD_TrueNegPionFromNeutralMeson_Pt",1000,0,25);
 659                         fTrueList[iCut]->Add(fHistoTrueNegPionFromNeutralMesonPt[iCut]);
 660                         fHistoTruePosPionFromNeutralMesonPt[iCut] = new TH1F("ESD_TruePosPionFromNeutralMeson_Pt","ESD_TruePosPionFromNeutralMeson_Pt",1000,0,25);
 661                         fTrueList[iCut]->Add(fHistoTruePosPionFromNeutralMesonPt[iCut]);
 662
 663                         fHistoTrueMotherPiPlPiMiPiZeroInvMassPt[iCut] = new TH2F("ESD_TrueMotherPiPlPiMiPiZero_InvMass_Pt","ESD_TrueMotherPiPlPiMiPiZero_InvMass_Pt",500,0.4,0.9,250,0,25);
 664                         fHistoTrueMotherPiPlPiMiPiZeroInvMassPt[iCut]->Sumw2();
 665                         fTrueList[iCut]->Add(fHistoTrueMotherPiPlPiMiPiZeroInvMassPt[iCut]);
 666
 667                         fHistoTrueMotherGammaGammaInvMassPt[iCut] = new TH2F("ESD_TrueMotherGG_InvMass_Pt","ESD_TrueMotherGG_InvMass_Pt",450,0.,0.45,250,0,25);
 668                         fTrueList[iCut]->Add(fHistoTrueMotherGammaGammaInvMassPt[iCut]);
 669                         fHistoTrueMotherGammaGammaFromEtaInvMassPt[iCut] = new TH2F("ESD_TrueMotherGGFromEta_InvMass_Pt","ESD_TrueMotherGGFromEta_InvMass_Pt",450,0.,0.45,250,0,25);
 670                         fTrueList[iCut]->Add(fHistoTrueMotherGammaGammaFromEtaInvMassPt[iCut]);
 671                         fHistoTrueMotherGammaGammaFromOmegaInvMassPt[iCut] = new TH2F("ESD_TrueMotherGGFromOmega_InvMass_Pt","ESD_TrueMotherGGFromOmega_InvMass_Pt",450,0.,0.45,250,0,25);
 672                         fTrueList[iCut]->Add(fHistoTrueMotherGammaGammaFromOmegaInvMassPt[iCut]);
 673
 674
 675 //                      if (fDoMesonQA){
 676                                 fHistoTruePionPionInvMassPt[iCut] = new TH2F("ESD_TruePiPlusPiNeg_InvMassPt","ESD_TruePiPlusPiNeg_InvMassPt",2000,0.,2.,200,0.,20.);
 677                                 fTrueList[iCut]->Add(fHistoTruePionPionInvMassPt[iCut]);
 678                                 fHistoTruePionPionFromSameMotherInvMassPt[iCut] = new TH2F("ESD_TruePiPlusPiNegFromSameMother_InvMassPt","ESD_TruePiPlusPiNegFromSameMother_InvMassPt",2000,0.,2.,200,0.,20.);
 679                                 fTrueList[iCut]->Add(fHistoTruePionPionFromSameMotherInvMassPt[iCut]);
 680                                 fHistoTruePionPionFromEtaInvMassPt[iCut] = new TH2F("ESD_TruePiPlusPiNegFromEta_InvMassPt","ESD_TruePiPlusPiNegFromEta_InvMassPt",2000,0.,2.,200,0.,20.);
 681                                 fTrueList[iCut]->Add(fHistoTruePionPionFromEtaInvMassPt[iCut]);
 682                                 fHistoTruePionPionFromOmegaInvMassPt[iCut] = new TH2F("ESD_TruePiPlusPiNegFromOmega_InvMassPt","ESD_TruePiPlusPiNegFromOmega_InvMassPt",2000,0.,2.,200,0.,20.);
 683                                 fTrueList[iCut]->Add(fHistoTruePionPionFromOmegaInvMassPt[iCut]);
 684
 685                                 //                      }
 686                 }
 687         }
 688
 689
 690
 691         InitBack(); // Init Background Handler
 692
 693         fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
 694         if(!fV0Reader){printf("Error: No V0 Reader");return;} // GetV0Reader
 695
 696         if(fV0Reader)
 697                 if((AliConversionPhotonCuts*)fV0Reader->GetConversionCuts())
 698                         if(((AliConversionPhotonCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms())
 699                                 fOutputContainer->Add(((AliConversionPhotonCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms());
 700
 701
 702
 703         fPionSelector=(AliPrimaryPionSelector*)AliAnalysisManager::GetAnalysisManager()->GetTask("PionSelector");
 704         if(!fPionSelector){printf("Error: No PionSelector");return;} // GetV0Reader
 705
 706         if( fPionSelector ){
 707                 if ( ((AliPrimaryPionCuts*)fPionSelector->GetPrimaryPionCuts())->GetCutHistograms() ){
 708                         fOutputContainer->Add( ((AliPrimaryPionCuts*)fPionSelector->GetPrimaryPionCuts())->GetCutHistograms() );
 709                 }
 710         }
 711
 712         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 713                 if( fPionCutArray ){
 714                         if( ((AliPrimaryPionCuts*)fPionCutArray->At(iCut))->GetCutHistograms() ) {
 715                                 fCutFolder[iCut]->Add( ((AliPrimaryPionCuts*)fPionCutArray->At(iCut))->GetCutHistograms() );
 716                         }
 717                 }
 718                 if (fNeutralPionMode < 2){
 719                         if( fGammaCutArray ) {
 720                                 if( ((AliConversionPhotonCuts*)fGammaCutArray->At(iCut))->GetCutHistograms() ) {
 721                                         fCutFolder[iCut]->Add( ((AliConversionPhotonCuts*)fGammaCutArray->At(iCut))->GetCutHistograms()  );
 722                                 }
 723                         }
 724                 }
 725                 if (fNeutralPionMode > 0){
 726                         if( fClusterCutArray ) {
 727                                 if( ((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutHistograms() ) {
 728                                         fCutFolder[iCut]->Add( ((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutHistograms()  );
 729                                 }
 730                         }
 731                 }
 732                 if( fNeutralPionMesonCutArray  ) {
 733                         if( ((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(iCut))->GetCutHistograms() ) {
 734                                 fCutFolder[iCut]->Add( ((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(iCut))->GetCutHistograms());
 735                         }
 736                 }
 737                 if( fMesonCutArray  ) {
 738                         if( ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutHistograms() ) {
 739                                 fCutFolder[iCut]->Add( ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutHistograms());
 740                         }
 741                 }
 742         }
 743
 744         PostData(1, fOutputContainer);
 745
 746 }
 747
 748 //______________________________________________________________________
 749 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::UserExec(Option_t *){
 750
 751         //
 752         // Execute analysis for current event
 753         //
 754
 755         fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
 756         if(!fV0Reader){printf("Error: No V0 Reader");return;} // GetV0Reader
 757
 758         Int_t eventQuality = ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEventQuality();
 759
 760         if(eventQuality == 2 || eventQuality == 3){// Event Not Accepted due to MC event missing or wrong trigger for V0ReaderV1
 761                 for(Int_t iCut = 0; iCut<fnCuts; iCut++){
 762                         fHistoNEvents[iCut]->Fill(eventQuality);
 763                 }
 764                 return;
 765         }
 766
 767         fPionSelector=(AliPrimaryPionSelector*)AliAnalysisManager::GetAnalysisManager()->GetTask("PionSelector");
 768         if(!fPionSelector){printf("Error: No PionSelector");return;} // GetV0Reader
 769
 770         if(fIsMC) fMCEvent     =  MCEvent();
 771         fESDEvent        = (AliESDEvent*)InputEvent();
 772         fReaderGammas    = fV0Reader->GetReconstructedGammas(); // Gammas from default Cut
 773         fSelectorNegPionIndex = fPionSelector->GetReconstructedNegPionIndex(); // Electrons from default Cut
 774         fSelectorPosPionIndex = fPionSelector->GetReconstructedPosPionIndex(); // Positrons from default Cut
 775
 776         fNumberOfESDTracks = fV0Reader->GetNumberOfPrimaryTracks();
 777         //AddTaskContainers(); //Add conatiner
 778
 779         for(Int_t iCut = 0; iCut<fnCuts; iCut++){
 780                 fiCut = iCut;
 781                 Int_t eventNotAccepted = ((AliConvEventCuts*)fEventCutArray->At(iCut))->IsEventAcceptedByCut(fV0Reader->GetEventCuts(),fInputEvent,fMCEvent,fIsHeavyIon);
 782
 783                 if(eventNotAccepted){
 784                         //                      cout << "event rejected due to wrong trigger: " <<eventNotAccepted << endl;
 785                         fHistoNEvents[iCut]->Fill(eventNotAccepted); // Check Centrality, PileUp, SDD and V0AND --> Not Accepted => eventQuality = 1
 786                         continue;
 787                 }
 788
 789                 if(eventQuality != 0){// Event Not Accepted
 790                         //                      cout << "event rejected due to: " <<eventQuality << endl;
 791                         fHistoNEvents[iCut]->Fill(eventQuality);
 792                         continue;
 793                 }
 794
 795                 fHistoNEvents[iCut]->Fill(eventQuality);
 796                 fHistoNGoodESDTracks[iCut]->Fill(fNumberOfESDTracks);
 797
 798                 if(fMCEvent){ // Process MC Particle
 799                         fMCStack = fMCEvent->Stack();
 800                         if(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetSignalRejection() != 0){
 801                                 ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetNotRejectedParticles(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetSignalRejection(),
 802                                                                                                                                                                                 ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetAcceptedHeader(),
 803                                                                                                                                                                                 fMCEvent);
 804                         }
 805                         ProcessMCParticles();
 806                 }
 807
 808                 if (fNeutralPionMode < 2){
 809                         ProcessConversionPhotonCandidates(); // Process this cuts conversion gammas
 810                 }
 811                 if (fNeutralPionMode > 0){
 812                         ProcessCaloPhotonCandidates(); // Process this cuts calo gammas
 813                 }
 814
 815                 if (fNeutralPionMode == 0 ){
 816                         ProcessNeutralPionCandidatesPureConversions(); // Process neutral pion candidates purely from conversions
 817                 }
 818                 if (fNeutralPionMode == 1){
 819                         ProcessNeutralPionCandidatesMixedConvCalo(); // Process neutral pion candidates mixed conv and calo
 820                 }
 821                 if (fNeutralPionMode == 2){
 822                         ProcessNeutralPionCandidatesPureCalo(); // Process neutral pion candidates purely from calo
 823                 }
 824
 825                 ProcessPionCandidates(); // Process this cuts gammas
 826
 827                 CalculateMesonCandidates();
 828                 CalculateBackground();
 829                 UpdateEventByEventData();
 830
 831
 832                 if (fGoodConvGammas->GetEntries()>0) fGoodConvGammas->Clear();
 833                 if (fClusterCandidates->GetEntries()>0) fClusterCandidates->Clear();
 834                 if (fNeutralPionCandidates->GetEntries()>0) fNeutralPionCandidates->Clear();
 835                 fGoodVirtualParticles->Clear(); // delete this cuts good gammas
 836         }
 837
 838         fSelectorNegPionIndex.clear();
 839         fSelectorPosPionIndex.clear();
 840
 841         PostData( 1, fOutputContainer );
 842 }
 843
 844 Bool_t AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::Notify(){
 845         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 846                 if( !((AliConvEventCuts*)fEventCutArray->At(iCut))->GetDoEtaShift() ){
 847                         fProfileEtaShift[iCut]->Fill(0.,0.);
 848                         continue; // No Eta Shift requested, continue
 849                 }
 850                 if( ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift() == 0.0){ // Eta Shift requested but not set, get shift automatically
 851                         ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCorrectEtaShiftFromPeriod(fV0Reader->GetPeriodName());
 852                         ((AliConvEventCuts*)fEventCutArray->At(iCut))->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
 853                         ((AliPrimaryPionCuts*)fPionCutArray->At(iCut))->SetEtaShift( ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift() );
 854                         fProfileEtaShift[iCut]->Fill(0.,(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift()));
 855                         continue;
 856                 } else {
 857                         printf(" Eta t PiPlusPiMinus Gamma Task %s :: Eta Shift Manually Set to %f \n\n",
 858                         (((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber()).Data(),((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift());
 859                         ((AliConvEventCuts*)fEventCutArray->At(iCut))->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
 860                         ((AliPrimaryPionCuts*)fPionCutArray->At(iCut))->SetEtaShift( ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift() );
 861                         fProfileEtaShift[iCut]->Fill(0.,(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift()));
 862                 }
 863         }
 864         return kTRUE;
 865 }
 866
 867
 868 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::Terminate(const Option_t *){
 869 ///Grid
 870 }
 871
 872
 873 //________________________________________________________________________
 874 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessCaloPhotonCandidates()
 875 {
 876
 877         Int_t nclus = 0;
 878         nclus = fInputEvent->GetNumberOfCaloClusters();
 879
 880 //      cout << nclus << endl;
 881
 882         if(nclus == 0)  return;
 883
 884         // vertex
 885         Double_t vertex[3] = {0};
 886         InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 887
 888         // Loop over EMCal clusters
 889         for(Long_t i = 0; i < nclus; i++){
 890
 891                 AliVCluster* clus = NULL;
 892                 clus = fInputEvent->GetCaloCluster(i);
 893                 if (!clus) continue;
 894                 if(!((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->ClusterIsSelected(clus,fInputEvent,fIsMC)) continue;
 895                 // TLorentzvector with cluster
 896                 TLorentzVector clusterVector;
 897                 clus->GetMomentum(clusterVector,vertex);
 898
 899                 TLorentzVector* tmpvec = new TLorentzVector();
 900                 tmpvec->SetPxPyPzE(clusterVector.Px(),clusterVector.Py(),clusterVector.Pz(),clusterVector.E());
 901
 902                 // convert to AODConversionPhoton
 903                 AliAODConversionPhoton *PhotonCandidate=new AliAODConversionPhoton(tmpvec);
 904                 if(!PhotonCandidate) continue;
 905
 906                 // Flag Photon as CaloPhoton
 907                 PhotonCandidate->SetIsCaloPhoton();
 908                 PhotonCandidate->SetCaloClusterRef(i);
 909                 // get MC label
 910                 if(fIsMC){
 911                         Int_t* mclabelsCluster = clus->GetLabels();
 912                         PhotonCandidate->SetNCaloPhotonMCLabels(clus->GetNLabels());
 913 //                      cout << clus->GetNLabels() << endl;
 914                         if (clus->GetNLabels()>0){
 915                                 for (Int_t k =0; k< (Int_t)clus->GetNLabels(); k++){
 916                                         if (k< 20)PhotonCandidate->SetCaloPhotonMCLabel(k,mclabelsCluster[k]);
 917 //                                      Int_t pdgCode = fMCStack->Particle(mclabelsCluster[k])->GetPdgCode();
 918 //                                      cout << "label " << k << "\t" << mclabelsCluster[k] << " pdg code: " << pdgCode << endl;
 919                                 }
 920                         }
 921                 }
 922
 923                 fIsFromMBHeader = kTRUE;
 924                 // test whether largest contribution to cluster orginates in added signals
 925                 if (fIsMC && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetCaloPhotonMCLabel(0), fMCStack, fInputEvent) == 0) fIsFromMBHeader = kFALSE;
 926
 927                 if (fIsFromMBHeader){
 928                         fHistoClusterGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
 929                         fHistoClusterGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
 930                 }
 931                 fClusterCandidates->Add(PhotonCandidate); // if no second loop is required add to events good gammas
 932
 933                 if(fIsMC){
 934 //                      if(fInputEvent->IsA()==AliESDEvent::Class()){
 935                                 ProcessTrueCaloPhotonCandidates(PhotonCandidate);
 936 //                      } else {
 937 //                              ProcessTrueClusterCandidatesAOD(PhotonCandidate);
 938 //                      }
 939                 }
 940
 941                 delete tmpvec;
 942         }
 943
 944 }
 945
 946 //________________________________________________________________________
 947 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessTrueCaloPhotonCandidates(AliAODConversionPhoton *TruePhotonCandidate)
 948 {
 949         TParticle *Photon = NULL;
 950         if (!TruePhotonCandidate->GetIsCaloPhoton()) AliFatal("CaloPhotonFlag has not been set task will abort");
 951 //      fHistoTrueNLabelsInClus[fiCut]->Fill(TruePhotonCandidate->GetNCaloPhotonMCLabels());
 952
 953         if (TruePhotonCandidate->GetNCaloPhotonMCLabels()>0)Photon = fMCStack->Particle(TruePhotonCandidate->GetCaloPhotonMCLabel(0));
 954                 else return;
 955
 956         if(Photon == NULL){
 957         //    cout << "no photon" << endl;
 958                 return;
 959         }
 960
 961 //      Int_t pdgCodeParticle = Photon->GetPdgCode();
 962         TruePhotonCandidate->SetCaloPhotonMCFlags(fMCStack);
 963
 964         // True Photon
 965         if(fIsFromMBHeader){
 966                 if(TruePhotonCandidate->GetCaloPhotonMCLabel(0) <= fMCStack->GetNprimary()){
 967                         if (TruePhotonCandidate->IsLargestComponentPhoton()){
 968                                 fHistoTrueClusterGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
 969                                 if (GammaIsNeutralMesonPiPlPiMiPiZeroDaughter(TruePhotonCandidate->GetCaloPhotonMCLabel(0))){
 970                                         fHistoTrueClusterGammaFromNeutralMesonPt[fiCut]->Fill(TruePhotonCandidate->Pt());
 971                                 }
 972                         }
 973                         if (TruePhotonCandidate->IsLargestComponentElectron() && TruePhotonCandidate->IsConversion()){
 974                                         fHistoTrueClusterGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
 975                                         if (GammaIsNeutralMesonPiPlPiMiPiZeroDaughter(TruePhotonCandidate->GetCaloPhotonMCLabel(0))){
 976                                         fHistoTrueClusterGammaFromNeutralMesonPt[fiCut]->Fill(TruePhotonCandidate->Pt());
 977                                 }
 978                         }
 979                 }
 980         }
 981         return;
 982 }
 983
 984
 985
 986 //________________________________________________________________________
 987 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessConversionPhotonCandidates(){
 988         Int_t nV0 = 0;
 989         TList *GoodGammasStepOne = new TList();
 990         TList *GoodGammasStepTwo = new TList();
 991         // Loop over Photon Candidates allocated by ReaderV1
 992
 993         for(Int_t i = 0; i < fReaderGammas->GetEntriesFast(); i++){
 994                 AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) fReaderGammas->At(i);
 995                 if(!PhotonCandidate) continue;
 996
 997                 fIsFromMBHeader = kTRUE;
 998
 999                 if( fMCEvent && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0 ){
1000                         Int_t isPosFromMBHeader
1001                                 = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack, fInputEvent);
1002                         if(isPosFromMBHeader == 0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1003                         Int_t isNegFromMBHeader
1004                                 = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack,fInputEvent);
1005                         if(isNegFromMBHeader == 0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1006                         if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE;
1007                 }
1008
1009                 if(!((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->PhotonIsSelected(PhotonCandidate,fESDEvent)) continue;
1010
1011                 if(!((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->UseElecSharingCut() &&
1012                         !((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->UseToCloseV0sCut()){ // if no post reader loop is required add to events good gammas
1013
1014                         fGoodConvGammas->Add(PhotonCandidate);
1015
1016                         if(fIsFromMBHeader){
1017                                 fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
1018                                 fHistoConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
1019                         }
1020
1021                         if(fMCEvent){
1022                                 ProcessTrueConversionPhotonCandidates(PhotonCandidate);
1023                         }
1024                 } else if(((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->UseElecSharingCut()){ // if Shared Electron cut is enabled, Fill array, add to step one
1025                         ((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->FillElectonLabelArray(PhotonCandidate,nV0);
1026                         nV0++;
1027                         GoodGammasStepOne->Add(PhotonCandidate);
1028                 } else if(!((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->UseElecSharingCut() &&
1029                                 ((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->UseToCloseV0sCut()){ // shared electron is disabled, step one not needed -> step two
1030                         GoodGammasStepTwo->Add(PhotonCandidate);
1031                 }
1032         }
1033
1034
1035         if(((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->UseElecSharingCut()){
1036                 for(Int_t i = 0;i<GoodGammasStepOne->GetEntries();i++){
1037                         AliAODConversionPhoton *PhotonCandidate= (AliAODConversionPhoton*) GoodGammasStepOne->At(i);
1038                         if(!PhotonCandidate) continue;
1039                         fIsFromMBHeader = kTRUE;
1040                         if(fMCEvent && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0){
1041                                 Int_t isPosFromMBHeader
1042                                 = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack,fInputEvent);
1043                                 Int_t isNegFromMBHeader
1044                                 = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack,fInputEvent);
1045                                 if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE;
1046                         }
1047                         if(!((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->RejectSharedElectronV0s(PhotonCandidate,i,GoodGammasStepOne->GetEntries())) continue;
1048                         if(!((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->UseToCloseV0sCut()){ // To Colse v0s cut diabled, step two not needed
1049                                 fGoodConvGammas->Add(PhotonCandidate);
1050                                 if(fIsFromMBHeader){
1051                                         fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
1052                                         fHistoConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
1053                                 }
1054                                 if(fMCEvent){
1055                                         ProcessTrueConversionPhotonCandidates(PhotonCandidate);
1056                                 }
1057                         }
1058                         else GoodGammasStepTwo->Add(PhotonCandidate); // Close v0s cut enabled -> add to list two
1059                 }
1060         }
1061         if(((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->UseToCloseV0sCut()){
1062                 for(Int_t i = 0;i<GoodGammasStepTwo->GetEntries();i++){
1063                         AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) GoodGammasStepTwo->At(i);
1064                         if(!PhotonCandidate) continue;
1065
1066                         if(fMCEvent && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0){
1067                                 Int_t isPosFromMBHeader
1068                                 = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack,fInputEvent);
1069                                 Int_t isNegFromMBHeader
1070                                 = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack,fInputEvent);
1071                                 if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE;
1072                         }
1073
1074                         if(!((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->RejectToCloseV0s(PhotonCandidate,GoodGammasStepTwo,i)) continue;
1075                         fGoodConvGammas->Add(PhotonCandidate); // Add gamma to current cut TList
1076
1077                         if(fIsFromMBHeader){
1078                                 fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt()); // Differences to old V0Reader in p_t due to conversion KF->TLorentzVector
1079                                 fHistoConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
1080                         }
1081
1082                         if(fMCEvent){
1083                                 ProcessTrueConversionPhotonCandidates(PhotonCandidate);
1084                         }
1085                 }
1086         }
1087
1088         delete GoodGammasStepOne;
1089         GoodGammasStepOne = 0x0;
1090         delete GoodGammasStepTwo;
1091         GoodGammasStepTwo = 0x0;
1092 }
1093
1094 //________________________________________________________________________
1095 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessTrueConversionPhotonCandidates(AliAODConversionPhoton *TruePhotonCandidate)
1096 {
1097         // Process True Photons
1098         AliStack *MCStack = fMCEvent->Stack();
1099         TParticle *posDaughter = TruePhotonCandidate->GetPositiveMCDaughter(MCStack);
1100         TParticle *negDaughter = TruePhotonCandidate->GetNegativeMCDaughter(MCStack);
1101
1102         if(posDaughter == NULL || negDaughter == NULL) return; // One particle does not exist
1103         if(posDaughter->GetMother(0) != negDaughter->GetMother(0)){  // Not Same Mother == Combinatorial Bck
1104                 return;
1105         }
1106
1107         else if (posDaughter->GetMother(0) == -1){
1108                 return;
1109         }
1110
1111         if(TMath::Abs(posDaughter->GetPdgCode())!=11 || TMath::Abs(negDaughter->GetPdgCode())!=11) return; //One Particle is not electron
1112         if(posDaughter->GetPdgCode()==negDaughter->GetPdgCode()) return; // Same Charge
1113         if(posDaughter->GetUniqueID() != 5 || negDaughter->GetUniqueID() !=5) return;// check if the daughters come from a conversion
1114
1115         TParticle *Photon = TruePhotonCandidate->GetMCParticle(MCStack);
1116         if(Photon->GetPdgCode() != 22) return; // Mother is no Photon
1117
1118         // True Photon
1119
1120         Int_t labelGamma = TruePhotonCandidate->GetMCParticleLabel(MCStack);
1121
1122         if( labelGamma < MCStack->GetNprimary() ){
1123                 if( fIsFromMBHeader ){
1124                         fHistoTrueConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1125                         if (GammaIsNeutralMesonPiPlPiMiPiZeroDaughter(labelGamma)){
1126                                 fHistoTrueConvGammaFromNeutralMesonPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1127                         }
1128                 }
1129         }
1130 }
1131
1132 //________________________________________________________________________
1133 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessNeutralPionCandidatesPureConversions(){
1134         // Conversion Gammas
1135         if(fGoodConvGammas->GetEntries()>1){
1136                 for(Int_t firstGammaIndex=0;firstGammaIndex<fGoodConvGammas->GetEntries()-1;firstGammaIndex++){
1137                         AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGoodConvGammas->At(firstGammaIndex));
1138                         if (gamma0==NULL) continue;
1139                         for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fGoodConvGammas->GetEntries();secondGammaIndex++){
1140                                 AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fGoodConvGammas->At(secondGammaIndex));
1141                                 //Check for same Electron ID
1142                                 if (gamma1==NULL) continue;
1143                                 if(gamma0->GetTrackLabelPositive() == gamma1->GetTrackLabelPositive() ||
1144                                 gamma0->GetTrackLabelNegative() == gamma1->GetTrackLabelNegative() ||
1145                                 gamma0->GetTrackLabelNegative() == gamma1->GetTrackLabelPositive() ||
1146                                 gamma0->GetTrackLabelPositive() == gamma1->GetTrackLabelNegative() ) continue;
1147
1148                                 AliAODConversionMother *pi0cand = new AliAODConversionMother(gamma0,gamma1);
1149                                 pi0cand->SetLabels(firstGammaIndex,secondGammaIndex);
1150
1151                                 pi0cand->CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
1152                                 if((((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(fiCut))->MesonIsSelected(pi0cand,kTRUE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift()))){
1153                                         fHistoGammaGammaInvMassPt[fiCut]->Fill(pi0cand->M(),pi0cand->Pt());
1154                                         if (pi0cand->M() > ((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(fiCut))->GetSelectionLow() && pi0cand->M() < ((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(fiCut))->GetSelectionHigh()){
1155                                                 fNeutralPionCandidates->Add(pi0cand);
1156 //                                              cout << "Pi0 candidate " << pi0cand->M() << "\t" << pi0cand->Pt() << endl;
1157                                         }
1158
1159                                         if(fIsMC){
1160                                                 if(fInputEvent->IsA()==AliESDEvent::Class())
1161                                                         ProcessTrueNeutralPionCandidatesPureConversions(pi0cand,gamma0,gamma1);
1162                                                 if(fInputEvent->IsA()==AliAODEvent::Class())
1163                                                         ProcessTrueNeutralPionCandidatesPureConversionsAOD(pi0cand,gamma0,gamma1);
1164                                         }
1165                                 }
1166 //                              delete pi0cand;
1167 //                              pi0cand=0x0;
1168                         }
1169                 }
1170         }
1171 }
1172
1173
1174 //________________________________________________________________________
1175 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessNeutralPionCandidatesPureCalo(){
1176
1177         // Conversion Gammas
1178         if(fClusterCandidates->GetEntries()>0){
1179
1180                 // vertex
1181                 Double_t vertex[3] = {0};
1182                 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1183
1184                 for(Int_t firstGammaIndex=0;firstGammaIndex<fClusterCandidates->GetEntries();firstGammaIndex++){
1185                         AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fClusterCandidates->At(firstGammaIndex));
1186                         if (gamma0==NULL) continue;
1187
1188                         for(Int_t secondGammaIndex=0;secondGammaIndex<fClusterCandidates->GetEntries();secondGammaIndex++){
1189                                 if (firstGammaIndex == secondGammaIndex) continue;
1190                                 AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fClusterCandidates->At(secondGammaIndex));
1191                                 if (gamma1==NULL) continue;
1192
1193                                 AliAODConversionMother *pi0cand = new AliAODConversionMother(gamma0,gamma1);
1194                                 pi0cand->SetLabels(firstGammaIndex,secondGammaIndex);
1195
1196                                 if((((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(fiCut))->MesonIsSelected(pi0cand,kTRUE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift()))){
1197                                         fHistoGammaGammaInvMassPt[fiCut]->Fill(pi0cand->M(),pi0cand->Pt());
1198                                         if (pi0cand->M() > ((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(fiCut))->GetSelectionLow() && pi0cand->M() < ((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(fiCut))->GetSelectionHigh()){
1199                                                 fNeutralPionCandidates->Add(pi0cand);
1200 //                                              cout << "Pi0 candidate " << pi0cand->M() << "\t" << pi0cand->Pt() << endl;
1201                                         }
1202
1203                                         if(fIsMC){
1204 //                                              if(fInputEvent->IsA()==AliESDEvent::Class())
1205                                                         ProcessTrueNeutralPionCandidatesPureCalo(pi0cand,gamma0,gamma1);
1206 //                                              if(fInputEvent->IsA()==AliAODEvent::Class())
1207 //                                                      ProcessTrueNeutralPionCandidatesPureConversionsAOD(pi0cand,gamma0,gamma1);
1208                                         }
1209                                 }
1210                         }
1211                 }
1212         }
1213 }
1214
1215 //______________________________________________________________________
1216 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessTrueNeutralPionCandidatesPureCalo( AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1)
1217 {
1218         // Process True Mesons
1219         AliStack *MCStack = fMCEvent->Stack();
1220
1221         Bool_t isTruePi0 = kFALSE;
1222         Int_t gamma0MCLabel = TrueGammaCandidate0->GetCaloPhotonMCLabel(0);     // get most probable MC label
1223         Int_t gamma0MotherLabel = -1;
1224         Int_t motherRealLabel = -1;
1225
1226         if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1227                 TParticle * gammaMC0 = (TParticle*)MCStack->Particle(gamma0MCLabel);
1228                 if (TrueGammaCandidate0->IsLargestComponentPhoton() || TrueGammaCandidate0->IsLargestComponentElectron()){              // largest component is electro magnetic
1229                         // get mother of interest (pi0 or eta)
1230                         if (TrueGammaCandidate0->IsLargestComponentPhoton()){                                                                                                           // for photons its the direct mother
1231                                 gamma0MotherLabel=gammaMC0->GetMother(0);
1232                                 motherRealLabel=gammaMC0->GetFirstMother();
1233                         } else if (TrueGammaCandidate0->IsLargestComponentElectron()){                                                                                          // for electrons its either the direct mother or for conversions the grandmother
1234                                 if (TrueGammaCandidate0->IsConversion()){
1235                                         gamma0MotherLabel=MCStack->Particle(gammaMC0->GetMother(0))->GetMother(0);
1236                                         motherRealLabel=MCStack->Particle(gammaMC0->GetMother(0))->GetMother(0);
1237                                 } else {
1238                                         gamma0MotherLabel=gammaMC0->GetMother(0);
1239                                         motherRealLabel=gammaMC0->GetMother(0);
1240                                 }
1241                         }
1242                 }
1243         }
1244
1245         if (!TrueGammaCandidate1->GetIsCaloPhoton()) AliFatal("CaloPhotonFlag has not been set. Aborting");
1246
1247         Int_t gamma1MCLabel = TrueGammaCandidate1->GetCaloPhotonMCLabel(0);     // get most probable MC label
1248         Int_t gamma1MotherLabel = -1;
1249         // check if
1250         if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1251                 // Daughters Gamma 1
1252                 TParticle * gammaMC1 = (TParticle*)MCStack->Particle(gamma1MCLabel);
1253                 if (TrueGammaCandidate1->IsLargestComponentPhoton() || TrueGammaCandidate1->IsLargestComponentElectron()){              // largest component is electro magnetic
1254                         // get mother of interest (pi0 or eta)
1255                         if (TrueGammaCandidate1->IsLargestComponentPhoton()){                                                                                                           // for photons its the direct mother
1256                                 gamma1MotherLabel=gammaMC1->GetMother(0);
1257                         } else if (TrueGammaCandidate1->IsLargestComponentElectron()){                                                                                          // for electrons its either the direct mother or for conversions the grandmother
1258                                 if (TrueGammaCandidate1->IsConversion()) gamma1MotherLabel=MCStack->Particle(gammaMC1->GetMother(0))->GetMother(0);
1259                                 else gamma1MotherLabel=gammaMC1->GetMother(0);
1260                         }
1261                 }
1262         }
1263
1264         if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
1265                 if(((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetPdgCode() == 111){
1266                         isTruePi0=kTRUE;
1267                 }
1268         }
1269
1270         if(isTruePi0){// True Pion
1271                 Pi0Candidate->SetTrueMesonValue(1);
1272                 Pi0Candidate->SetMCLabel(motherRealLabel);
1273                 fHistoTrueMotherGammaGammaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1274                 if( IsEtaPiPlPiMiPiZeroDaughter(motherRealLabel) ) {
1275                         fHistoTrueMotherGammaGammaFromEtaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1276                 }
1277                 if( IsOmegaPiPlPiMiPiZeroDaughter(motherRealLabel) ) {
1278                         fHistoTrueMotherGammaGammaFromOmegaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1279                 }
1280         }
1281 }
1282
1283
1284
1285 //______________________________________________________________________
1286 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessTrueNeutralPionCandidatesPureConversions(AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1)
1287 {
1288         // Process True Mesons
1289         AliStack *MCStack = fMCEvent->Stack();
1290         if(TrueGammaCandidate0->GetV0Index()<fInputEvent->GetNumberOfV0s()){
1291                 Bool_t isTruePi0 = kFALSE;
1292                 Bool_t isTruePi0Dalitz = kFALSE;
1293                 Bool_t gamma0DalitzCand = kFALSE;
1294                 Bool_t gamma1DalitzCand = kFALSE;
1295                 Int_t gamma0MCLabel = TrueGammaCandidate0->GetMCParticleLabel(MCStack);
1296                 Int_t gamma0MotherLabel = -1;
1297                 Int_t motherRealLabel = -1;
1298                 if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1299                         // Daughters Gamma 0
1300                         TParticle * negativeMC = (TParticle*)TrueGammaCandidate0->GetNegativeMCDaughter(MCStack);
1301                         TParticle * positiveMC = (TParticle*)TrueGammaCandidate0->GetPositiveMCDaughter(MCStack);
1302                         TParticle * gammaMC0 = (TParticle*)MCStack->Particle(gamma0MCLabel);
1303                         if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1304                                 if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
1305                                         if(gammaMC0->GetPdgCode() == 22){ // ... with Gamma Mother
1306                                                 gamma0MotherLabel=gammaMC0->GetFirstMother();
1307                                                 motherRealLabel=gammaMC0->GetFirstMother();
1308                                         }
1309                                 }
1310                                 if(gammaMC0->GetPdgCode() ==111){ // Dalitz candidate
1311                                         gamma0DalitzCand = kTRUE;
1312                                         gamma0MotherLabel=-111;
1313                                         motherRealLabel=gamma0MCLabel;
1314                                 }
1315                         }
1316                 }
1317                 if(TrueGammaCandidate1->GetV0Index()<fInputEvent->GetNumberOfV0s()){
1318                         Int_t gamma1MCLabel = TrueGammaCandidate1->GetMCParticleLabel(MCStack);
1319                         Int_t gamma1MotherLabel = -1;
1320                         if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1321                                 // Daughters Gamma 1
1322                                 TParticle * negativeMC = (TParticle*)TrueGammaCandidate1->GetNegativeMCDaughter(MCStack);
1323                                 TParticle * positiveMC = (TParticle*)TrueGammaCandidate1->GetPositiveMCDaughter(MCStack);
1324                                 TParticle * gammaMC1 = (TParticle*)MCStack->Particle(gamma1MCLabel);
1325                                 if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1326                                         if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
1327                                                 if(gammaMC1->GetPdgCode() == 22){ // ... with Gamma Mother
1328                                                         gamma1MotherLabel=gammaMC1->GetFirstMother();
1329                                                 }
1330                                         }
1331                                         if(gammaMC1->GetPdgCode() ==111 ){ // Dalitz candidate
1332                                                 gamma1DalitzCand = kTRUE;
1333                                                 gamma1MotherLabel=-111;
1334                                         }
1335                                 }
1336                         }
1337                         if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
1338                                 if(((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetPdgCode() == 111){
1339                                         isTruePi0=kTRUE;
1340                                 }
1341                         }
1342
1343                         //Identify Dalitz candidate
1344                         if (gamma1DalitzCand || gamma0DalitzCand){
1345                                 if (gamma0DalitzCand && gamma0MCLabel >=0 && gamma0MCLabel==gamma1MotherLabel){
1346                                         if (gamma0MotherLabel == -111) isTruePi0Dalitz = kTRUE;
1347                                 }
1348                                 if (gamma1DalitzCand && gamma1MCLabel >=0 && gamma1MCLabel==gamma0MotherLabel){
1349                                         if (gamma1MotherLabel == -111) isTruePi0Dalitz = kTRUE;
1350                                 }
1351                         }
1352
1353
1354                         if(isTruePi0 || isTruePi0Dalitz){// True Pion
1355                                 Pi0Candidate->SetTrueMesonValue(1);
1356                                 Pi0Candidate->SetMCLabel(motherRealLabel);
1357                                 fHistoTrueMotherGammaGammaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1358                                 if( IsEtaPiPlPiMiPiZeroDaughter(motherRealLabel) ) {
1359                                         fHistoTrueMotherGammaGammaFromEtaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1360                                 }
1361                                 if( IsOmegaPiPlPiMiPiZeroDaughter(motherRealLabel) ) {
1362                                         fHistoTrueMotherGammaGammaFromOmegaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1363                                 }
1364                         }
1365                 }
1366         }
1367 }
1368
1369 //______________________________________________________________________
1370 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessTrueNeutralPionCandidatesPureConversionsAOD(AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1)
1371 {
1372
1373         // Process True Mesons
1374         TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1375         Bool_t isTruePi0 = kFALSE;
1376         Bool_t isTruePi0Dalitz = kFALSE;
1377         Bool_t gamma0DalitzCand = kFALSE;
1378         Bool_t gamma1DalitzCand = kFALSE;
1379         Int_t motherRealLabel = -1;
1380
1381         if (AODMCTrackArray!=NULL && TrueGammaCandidate0 != NULL){
1382                 AliAODMCParticle *positiveMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate0->GetMCLabelPositive()));
1383                 AliAODMCParticle *negativeMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate0->GetMCLabelNegative()));
1384
1385                 Int_t gamma0MCLabel = -1;
1386                 Int_t gamma0MotherLabel = -1;
1387                 if(!positiveMC||!negativeMC)
1388                         return;
1389
1390                 if(positiveMC->GetMother()>-1&&(negativeMC->GetMother() == positiveMC->GetMother())){
1391                         gamma0MCLabel = positiveMC->GetMother();
1392                 }
1393
1394                 if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1395                         // Daughters Gamma 0
1396                         AliAODMCParticle * gammaMC0 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma0MCLabel));
1397                         if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1398                                 if(((positiveMC->GetMCProcessCode())) == 5 && ((negativeMC->GetMCProcessCode())) == 5){ // ... From Conversion ...
1399                                         if(gammaMC0->GetPdgCode() == 22){ // ... with Gamma Mother
1400                                                 gamma0MotherLabel=gammaMC0->GetMother();
1401                                                 motherRealLabel=gammaMC0->GetMother();
1402                                         }
1403                                 }
1404                                 if(gammaMC0->GetPdgCode() ==111){ // Dalitz candidate
1405                                         gamma0DalitzCand = kTRUE;
1406                                         gamma0MotherLabel=-111;
1407                                         motherRealLabel=gamma0MCLabel;
1408                                 }
1409                         }
1410                 }
1411                 positiveMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate1->GetMCLabelPositive()));
1412                 negativeMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate1->GetMCLabelNegative()));
1413
1414                 Int_t gamma1MCLabel = -1;
1415                 Int_t gamma1MotherLabel = -1;
1416                 if(!positiveMC||!negativeMC)
1417                         return;
1418
1419                 if(positiveMC->GetMother()>-1&&(negativeMC->GetMother() == positiveMC->GetMother())){
1420                         gamma1MCLabel = positiveMC->GetMother();
1421                 }
1422                 if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1423                         // Daughters Gamma 1
1424                         AliAODMCParticle * gammaMC1 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MCLabel));
1425                         if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1426                                 if(((positiveMC->GetMCProcessCode())) == 5 && ((negativeMC->GetMCProcessCode())) == 5){ // ... From Conversion ...
1427                                         if(gammaMC1->GetPdgCode() == 22){ // ... with Gamma Mother
1428                                         gamma1MotherLabel=gammaMC1->GetMother();
1429                                         }
1430                                 }
1431                                 if(gammaMC1->GetPdgCode() ==111 ){ // Dalitz candidate
1432                                                 gamma1DalitzCand = kTRUE;
1433                                                 gamma1MotherLabel=-111;
1434                                 }
1435                         }
1436                 }
1437                 if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
1438                         if(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->GetPdgCode() == 111){
1439                                 isTruePi0=kTRUE;
1440                         }
1441                 }
1442
1443                 //Identify Dalitz candidate
1444                 if (gamma1DalitzCand || gamma0DalitzCand){
1445                         if (gamma0DalitzCand && gamma0MCLabel >=0 && gamma0MCLabel==gamma1MotherLabel){
1446                                 if (gamma0MotherLabel == -111) isTruePi0Dalitz = kTRUE;
1447                         }
1448                         if (gamma1DalitzCand && gamma1MCLabel >=0 && gamma1MCLabel==gamma0MotherLabel){
1449                                 if (gamma1MotherLabel == -111) isTruePi0Dalitz = kTRUE;
1450                         }
1451                 }
1452
1453                 if(isTruePi0 || isTruePi0Dalitz){// True Pion
1454                         Pi0Candidate->SetTrueMesonValue(1);
1455                         Pi0Candidate->SetMCLabel(motherRealLabel);
1456                         fHistoTrueMotherGammaGammaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1457                         if( IsEtaPiPlPiMiPiZeroDaughter(motherRealLabel) ) {
1458                                 fHistoTrueMotherGammaGammaFromEtaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1459                         }
1460                         if( IsOmegaPiPlPiMiPiZeroDaughter(motherRealLabel) ) {
1461                                 fHistoTrueMotherGammaGammaFromOmegaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1462                         }
1463                 }
1464         }
1465         return;
1466 }
1467
1468
1469 //________________________________________________________________________
1470 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessNeutralPionCandidatesMixedConvCalo(){
1471
1472         // Conversion Gammas
1473         if(fGoodConvGammas->GetEntries()>0){
1474                 // vertex
1475                 Double_t vertex[3] = {0};
1476                 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1477
1478                 for(Int_t firstGammaIndex=0;firstGammaIndex<fGoodConvGammas->GetEntries();firstGammaIndex++){
1479                         AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGoodConvGammas->At(firstGammaIndex));
1480                         if (gamma0==NULL) continue;
1481
1482                         for(Int_t secondGammaIndex=0;secondGammaIndex<fClusterCandidates->GetEntries();secondGammaIndex++){
1483                                 Bool_t matched = kFALSE;
1484                                 AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fClusterCandidates->At(secondGammaIndex));
1485                                 if (gamma1==NULL) continue;
1486
1487                                 if (gamma1->GetIsCaloPhoton()){
1488                                         AliVCluster* cluster = fInputEvent->GetCaloCluster(gamma1->GetCaloClusterRef());
1489                                         matched = ((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->MatchConvPhotonToCluster(gamma0,cluster, fInputEvent );
1490                                 }
1491
1492                                 AliAODConversionMother *pi0cand = new AliAODConversionMother(gamma0,gamma1);
1493                                 pi0cand->SetLabels(firstGammaIndex,secondGammaIndex);
1494
1495                                 if((((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(fiCut))->MesonIsSelected(pi0cand,kTRUE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift()))){
1496                                         if (!matched){
1497                                                 fHistoGammaGammaInvMassPt[fiCut]->Fill(pi0cand->M(),pi0cand->Pt());
1498                                                 if (pi0cand->M() > ((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(fiCut))->GetSelectionLow() && pi0cand->M() < ((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(fiCut))->GetSelectionHigh()){
1499                                                         fNeutralPionCandidates->Add(pi0cand);
1500         //                                              cout << "Pi0 candidate " << pi0cand->M() << "\t" << pi0cand->Pt() << endl;
1501                                                 }
1502
1503                                                 if(fIsMC){
1504 //                                                      if(fInputEvent->IsA()==AliESDEvent::Class())
1505                                                                 ProcessTrueNeutralPionCandidatesMixedConvCalo(pi0cand,gamma0,gamma1);
1506 //                                                      if(fInputEvent->IsA()==AliAODEvent::Class())
1507 //                                                              ProcessTrueMesonCandidatesAOD(pi0cand,gamma0,gamma1, matched);
1508                                                 }
1509                                         }
1510                                 }
1511                         }
1512                 }
1513         }
1514 }
1515
1516 //______________________________________________________________________
1517 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessTrueNeutralPionCandidatesMixedConvCalo( AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1)
1518 {
1519         // Process True Mesons
1520         AliStack *MCStack = fMCEvent->Stack();
1521         if(TrueGammaCandidate0->GetV0Index()<fInputEvent->GetNumberOfV0s()){
1522                 Bool_t isTruePi0 = kFALSE;
1523                 Bool_t isTruePi0Dalitz = kFALSE;
1524                 Bool_t gamma0DalitzCand = kFALSE;
1525
1526                 Int_t gamma0MCLabel = TrueGammaCandidate0->GetMCParticleLabel(MCStack);
1527                 Int_t gamma0MotherLabel = -1;
1528                 Int_t motherRealLabel = -1;
1529                 if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1530                         // Daughters Gamma 0
1531                         TParticle * negativeMC = (TParticle*)TrueGammaCandidate0->GetNegativeMCDaughter(MCStack);
1532                         TParticle * positiveMC = (TParticle*)TrueGammaCandidate0->GetPositiveMCDaughter(MCStack);
1533                         TParticle * gammaMC0 = (TParticle*)MCStack->Particle(gamma0MCLabel);
1534                         if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1535                                 if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
1536                                         if(gammaMC0->GetPdgCode() == 22){ // ... with Gamma Mother
1537                                                 gamma0MotherLabel=gammaMC0->GetFirstMother();
1538                                                 motherRealLabel=gammaMC0->GetFirstMother();
1539                                         }
1540                                 }
1541                                 if(gammaMC0->GetPdgCode() ==111){ // Dalitz candidate
1542                                         gamma0DalitzCand = kTRUE;
1543                                         gamma0MotherLabel=-111;
1544                                         motherRealLabel=gamma0MCLabel;
1545                                 }
1546
1547                         }
1548                 }
1549
1550                 if (!TrueGammaCandidate1->GetIsCaloPhoton()) AliFatal("CaloPhotonFlag has not been set. Aborting");
1551
1552                 Int_t gamma1MCLabel = TrueGammaCandidate1->GetCaloPhotonMCLabel(0);     // get most probable MC label
1553                 Int_t gamma1MotherLabel = -1;
1554                 // check if
1555
1556                 if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1557                         // Daughters Gamma 1
1558                         TParticle * gammaMC1 = (TParticle*)MCStack->Particle(gamma1MCLabel);
1559                         if (TrueGammaCandidate1->IsLargestComponentPhoton() || TrueGammaCandidate1->IsLargestComponentElectron()){              // largest component is electro magnetic
1560                                 // get mother of interest (pi0 or eta)
1561                                 if (TrueGammaCandidate1->IsLargestComponentPhoton()){                                                                                                           // for photons its the direct mother
1562                                         gamma1MotherLabel=gammaMC1->GetMother(0);
1563                                 } else if (TrueGammaCandidate1->IsLargestComponentElectron()){                                                                                          // for electrons its either the direct mother or for conversions the grandmother
1564                                         if (TrueGammaCandidate1->IsConversion()) gamma1MotherLabel=MCStack->Particle(gammaMC1->GetMother(0))->GetMother(0);
1565                                         else gamma1MotherLabel=gammaMC1->GetMother(0);
1566                                 }
1567                         }
1568                 }
1569
1570                 if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
1571                         if(((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetPdgCode() == 111){
1572                                 isTruePi0=kTRUE;
1573                         }
1574                 }
1575
1576                 if (gamma0DalitzCand ){
1577                         if (gamma0DalitzCand && gamma0MCLabel >=0 && gamma0MCLabel==gamma1MotherLabel){
1578                                 if (gamma0MotherLabel == -111) isTruePi0Dalitz = kTRUE;
1579                         }
1580                 }
1581
1582                 if(isTruePi0 || isTruePi0Dalitz ){
1583                         Pi0Candidate->SetTrueMesonValue(1);
1584                         Pi0Candidate->SetMCLabel(motherRealLabel);
1585                         fHistoTrueMotherGammaGammaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1586                         if( IsEtaPiPlPiMiPiZeroDaughter(motherRealLabel) ) {
1587                                 fHistoTrueMotherGammaGammaFromEtaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1588                         }
1589                         if( IsOmegaPiPlPiMiPiZeroDaughter(motherRealLabel) ) {
1590                                 fHistoTrueMotherGammaGammaFromOmegaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1591                         }
1592                 }
1593         }
1594 }
1595
1596
1597
1598 //________________________________________________________________________
1599 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessPionCandidates(){
1600
1601         Double_t magField = fInputEvent->GetMagneticField();
1602         if( magField  < 0.0 ){
1603                 magField =  1.0;
1604         } else {
1605                 magField =  -1.0;
1606         }
1607
1608         vector<Int_t> lGoodNegPionIndexPrev(0);
1609         vector<Int_t> lGoodPosPionIndexPrev(0);
1610
1611         for(UInt_t i = 0; i < fSelectorNegPionIndex.size(); i++){
1612                 AliESDtrack* negPionCandidate = fESDEvent->GetTrack(fSelectorNegPionIndex[i]);
1613                 if(! ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->PionIsSelected(negPionCandidate) ) continue;
1614                 lGoodNegPionIndexPrev.push_back(   fSelectorNegPionIndex[i] );
1615                 fHistoNegPionPt[fiCut]->Fill(negPionCandidate->Pt());
1616                 fHistoNegPionPhi[fiCut]->Fill(negPionCandidate->Phi());
1617                 if( fMCEvent ) {
1618                         Int_t labelNegPion = TMath::Abs( negPionCandidate->GetLabel() );
1619                         if( labelNegPion < fMCStack->GetNtrack() ){
1620                                 TParticle* negPion = fMCStack->Particle(labelNegPion);
1621                                 if( negPion->GetPdgCode() ==  -211 ){
1622                                         if( labelNegPion < fMCStack->GetNprimary() ){
1623                                                 fHistoTrueNegPionPt[fiCut]->Fill(negPionCandidate->Pt());    //primary negPion
1624                                         }
1625                                         if( IsEtaPiPlPiMiPiZeroDaughter(labelNegPion) || IsOmegaPiPlPiMiPiZeroDaughter(labelNegPion) ) {
1626                                                 if( labelNegPion < fMCStack->GetNprimary() ) {
1627                                                         fHistoTrueNegPionFromNeutralMesonPt[fiCut]->Fill(negPionCandidate->Pt());
1628                                                 }
1629                                         }
1630                                 }
1631                         }
1632                 }
1633         }
1634
1635         for(UInt_t i = 0; i < fSelectorPosPionIndex.size(); i++){
1636                 AliESDtrack* posPionCandidate = fESDEvent->GetTrack( fSelectorPosPionIndex[i] );
1637                 if(! ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->PionIsSelected(posPionCandidate) ) continue;
1638                 lGoodPosPionIndexPrev.push_back(   fSelectorPosPionIndex[i]  );
1639                 fHistoPosPionPt[fiCut]->Fill( posPionCandidate->Pt() );
1640                 fHistoPosPionPhi[fiCut]->Fill( posPionCandidate->Phi() );
1641
1642                 if( fMCEvent ) {
1643                         Int_t labelPosPion = TMath::Abs( posPionCandidate->GetLabel() );
1644                         if( labelPosPion < fMCStack->GetNtrack() ) {
1645                                 TParticle* posPion = fMCStack->Particle(labelPosPion);
1646                                 if( posPion->GetPdgCode() ==  211 ){
1647                                         if( labelPosPion < fMCStack->GetNprimary() ){
1648                                                 fHistoTruePosPionPt[fiCut]->Fill(posPionCandidate->Pt());
1649                                         }
1650                                         if( IsEtaPiPlPiMiPiZeroDaughter(labelPosPion) || IsOmegaPiPlPiMiPiZeroDaughter(labelPosPion) ) {
1651                                                 if( labelPosPion < fMCStack->GetNprimary() ){
1652                                                         fHistoTruePosPionFromNeutralMesonPt[fiCut]->Fill(posPionCandidate->Pt());
1653                                                 }
1654                                         }
1655                                 }
1656                         }
1657                 }
1658         }
1659
1660
1661         for(UInt_t i = 0; i < lGoodNegPionIndexPrev.size(); i++){
1662
1663                 AliESDtrack *negPionCandidate = fESDEvent->GetTrack(lGoodNegPionIndexPrev[i]);
1664                 AliKFParticle negPionCandidateKF( *negPionCandidate->GetConstrainedParam(), 211 );
1665
1666                 for(UInt_t j = 0; j < lGoodPosPionIndexPrev.size(); j++){
1667
1668                         AliESDtrack *posPionCandidate = fESDEvent->GetTrack(lGoodPosPionIndexPrev[j]);
1669                         AliKFParticle posPionCandidateKF( *posPionCandidate->GetConstrainedParam(), 211 );
1670
1671                         AliKFConversionPhoton* virtualPhoton = NULL;
1672                         virtualPhoton = new AliKFConversionPhoton(negPionCandidateKF,posPionCandidateKF);
1673                         AliKFVertex primaryVertexImproved(*fInputEvent->GetPrimaryVertex());
1674 //                      primaryVertexImproved+=*virtualPhoton;
1675                         virtualPhoton->SetProductionVertex(primaryVertexImproved);
1676                         virtualPhoton->SetTrackLabels( lGoodPosPionIndexPrev[j], lGoodNegPionIndexPrev[i]);
1677
1678                         Int_t labeln=0;
1679                         Int_t labelp=0;
1680                         Int_t motherlabelp = 0;
1681                         Int_t motherlabeln = 0;
1682                         TParticle *fNegativeMCParticle =NULL;
1683                         TParticle *fPositiveMCParticle =NULL;
1684                         if( fMCEvent ) {
1685                                 labeln=TMath::Abs(negPionCandidate->GetLabel());
1686                                 labelp=TMath::Abs(posPionCandidate->GetLabel());
1687                                 fNegativeMCParticle = fMCStack->Particle(labeln);
1688                                 fPositiveMCParticle = fMCStack->Particle(labelp);
1689                                 motherlabeln = fNegativeMCParticle->GetMother(0);
1690                                 motherlabelp = fPositiveMCParticle->GetMother(0);
1691                                 if( fPositiveMCParticle && fNegativeMCParticle) {
1692                                         virtualPhoton->SetMCLabelPositive(labelp);
1693                                         virtualPhoton->SetMCLabelNegative(labeln);
1694                                 }
1695                         }
1696
1697                         AliAODConversionPhoton *vParticle = new AliAODConversionPhoton(virtualPhoton); //To apply mass 2 pion mass cut
1698                         if (((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->DoMassCut()){
1699                                 if (vParticle->GetMass() < ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->GetMassCut()){
1700                                         fGoodVirtualParticles->Add(  vParticle );
1701                                         fHistoPionPionInvMassPt[fiCut]->Fill( vParticle->GetMass(),vParticle->Pt());
1702                                 }
1703                         } else {
1704                                 fGoodVirtualParticles->Add(  vParticle );
1705                                 fHistoPionPionInvMassPt[fiCut]->Fill( vParticle->GetMass(),vParticle->Pt());
1706                         }
1707
1708                         Double_t clsToFPos = -1.0;
1709                         Double_t clsToFNeg = -1.0;
1710
1711                         Float_t dcaToVertexXYPos = -1.0;
1712                         Float_t dcaToVertexZPos  = -1.0;
1713                         Float_t dcaToVertexXYNeg = -1.0;
1714                         Float_t dcaToVertexZNeg  = -1.0;
1715
1716
1717
1718                         if ( fDoMesonQA ) {
1719                                 clsToFPos = ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->GetNFindableClustersTPC(posPionCandidate);
1720                                 clsToFNeg = ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->GetNFindableClustersTPC(negPionCandidate);
1721
1722                                 Float_t bPos[2];
1723                                 Float_t bCovPos[3];
1724                                 posPionCandidate->GetImpactParameters(bPos,bCovPos);
1725                                 if (bCovPos[0]<=0 || bCovPos[2]<=0) {
1726                                         AliDebug(1, "Estimated b resolution lower or equal zero!");
1727                                         bCovPos[0]=0; bCovPos[2]=0;
1728                                 }
1729
1730                                 Float_t bNeg[2];
1731                                 Float_t bCovNeg[3];
1732                                 posPionCandidate->GetImpactParameters(bNeg,bCovNeg);
1733                                 if (bCovNeg[0]<=0 || bCovNeg[2]<=0) {
1734                                         AliDebug(1, "Estimated b resolution lower or equal zero!");
1735                                         bCovNeg[0]=0; bCovNeg[2]=0;
1736                                 }
1737
1738                                 dcaToVertexXYPos = bPos[0];
1739                                 dcaToVertexZPos  = bPos[1];
1740                                 dcaToVertexXYNeg = bNeg[0];
1741                                 dcaToVertexZNeg  = bNeg[1];
1742
1743                                 fHistoNegPionEta[fiCut]->Fill( negPionCandidate->Eta() );
1744                                 fHistoPosPionEta[fiCut]->Fill( posPionCandidate->Eta() );
1745
1746                                 fHistoNegPionClsTPC[fiCut]->Fill(clsToFNeg,negPionCandidate->Pt());
1747                                 fHistoPosPionClsTPC[fiCut]->Fill(clsToFPos,posPionCandidate->Pt());
1748
1749                                 fHistoPionDCAxy[fiCut]->Fill(  dcaToVertexXYNeg, negPionCandidate->Pt() );
1750                                 fHistoPionDCAz[fiCut]->Fill(   dcaToVertexZNeg,  negPionCandidate->Pt() );
1751                                 fHistoPionDCAxy[fiCut]->Fill(  dcaToVertexXYPos, posPionCandidate->Pt() );
1752                                 fHistoPionDCAz[fiCut]->Fill(   dcaToVertexZPos,  posPionCandidate->Pt() );
1753
1754                                 fHistoPionTPCdEdxNSigma[fiCut]->Fill( posPionCandidate->P(),((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->GetPIDResponse()->NumberOfSigmasTPC(posPionCandidate, AliPID::kPion) );
1755                                 fHistoPionTPCdEdxNSigma[fiCut]->Fill( negPionCandidate->P(),((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->GetPIDResponse()->NumberOfSigmasTPC(negPionCandidate, AliPID::kPion) );
1756
1757                                 fHistoPionTPCdEdx[fiCut]->Fill( posPionCandidate->P(), TMath::Abs(posPionCandidate->GetTPCsignal()));
1758                                 fHistoPionTPCdEdx[fiCut]->Fill( negPionCandidate->P(), TMath::Abs(negPionCandidate->GetTPCsignal()));
1759                         }
1760
1761                         if (fMCEvent){
1762                                 if (fPositiveMCParticle && fNegativeMCParticle ) {
1763                                         if (((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->DoMassCut()){
1764                                                 if (vParticle->GetMass() < ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->GetMassCut()){
1765                                                         if(TMath::Abs(fNegativeMCParticle->GetPdgCode())==211 && TMath::Abs(fPositiveMCParticle->GetPdgCode())==211){  // Pions ...
1766                                                                 fHistoTruePionPionInvMassPt[fiCut]->Fill(vParticle->GetMass(),vParticle->Pt());
1767                                                                 if (motherlabeln == motherlabelp){
1768                                                                         fHistoTruePionPionFromSameMotherInvMassPt[fiCut]->Fill(vParticle->GetMass(),vParticle->Pt());
1769                                                                         if( IsEtaPiPlPiMiPiZeroDaughter(labeln) ) { //|| IsOmegaPiPlPiMiPiZeroDaughter(labeln)
1770                                                                                 fHistoTruePionPionFromEtaInvMassPt[fiCut]->Fill(vParticle->GetMass(),vParticle->Pt());
1771                                                                         }
1772                                                                         if( IsOmegaPiPlPiMiPiZeroDaughter(labeln) ) { //||
1773                                                                                 fHistoTruePionPionFromOmegaInvMassPt[fiCut]->Fill(vParticle->GetMass(),vParticle->Pt());
1774                                                                         }
1775                                                                 }
1776                                                         }
1777                                                 }
1778                                         } else {
1779                                                 if(TMath::Abs(fNegativeMCParticle->GetPdgCode())==211 && TMath::Abs(fPositiveMCParticle->GetPdgCode())==211){  // Pions ...
1780                                                         fHistoTruePionPionInvMassPt[fiCut]->Fill(vParticle->GetMass(),vParticle->Pt());
1781                                                         if (motherlabeln == motherlabelp){
1782                                                                 fHistoTruePionPionFromSameMotherInvMassPt[fiCut]->Fill(vParticle->GetMass(),vParticle->Pt());
1783                                                                 if( IsEtaPiPlPiMiPiZeroDaughter(labeln) ) { //|| IsOmegaPiPlPiMiPiZeroDaughter(labeln)
1784                                                                         fHistoTruePionPionFromEtaInvMassPt[fiCut]->Fill(vParticle->GetMass(),vParticle->Pt());
1785                                                                 }
1786                                                                 if( IsOmegaPiPlPiMiPiZeroDaughter(labeln) ) { //||
1787                                                                         fHistoTruePionPionFromOmegaInvMassPt[fiCut]->Fill(vParticle->GetMass(),vParticle->Pt());
1788                                                                 }
1789                                                         }
1790                                                 }
1791                                         }
1792                                 }
1793                         }
1794
1795                         delete virtualPhoton;
1796                         virtualPhoton=NULL;
1797
1798                 }
1799         }
1800 }
1801
1802 //_____________________________________________________________________________
1803 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessMCParticles(){
1804
1805         // Loop over all primary MC particle
1806
1807         for(Int_t i = 0; i < fMCStack->GetNprimary(); i++) {
1808
1809                 TParticle* particle = (TParticle *)fMCStack->Particle(i);
1810                 if (!particle) continue;
1811
1812                 Int_t isMCFromMBHeader = -1;
1813                 if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0){
1814                         isMCFromMBHeader
1815                                 = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent);
1816                         if(isMCFromMBHeader == 0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1817                 }
1818
1819                 if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack,fInputEvent)){
1820                         // find MC photons
1821                         if (fNeutralPionMode < 2){
1822                                 if(((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->PhotonIsSelectedMC(particle,fMCStack,kFALSE)){
1823                                         fHistoMCAllGammaPt[fiCut]->Fill(particle->Pt()); // All MC Gamma
1824                                         if(particle->GetMother(0) >-1){
1825                                                 if (fMCStack->Particle(particle->GetMother(0))->GetPdgCode() ==111){
1826                                                         if (fMCStack->Particle(particle->GetMother(0))->GetMother(0) > -1){
1827                                                                 if ( fMCStack->Particle((fMCStack->Particle(particle->GetMother(0)))->GetMother(0))->GetPdgCode() == 221 ||
1828                                                                         fMCStack->Particle((fMCStack->Particle(particle->GetMother(0)))->GetMother(0))->GetPdgCode() == 223 ){
1829                                                                         if ( fMCStack->Particle(particle->GetMother(0))->GetNDaughters()==3 )
1830                                                                                 fHistoMCGammaFromNeutralMesonPt[fiCut]->Fill(particle->Pt()); // All photons from eta or omega via pi0
1831                                                                 }
1832                                                         }
1833                                                 }
1834                                         }
1835                                 }
1836                         } else if (fNeutralPionMode == 2){
1837                                 if(((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->ClusterIsSelectedMC(particle,fMCStack)){
1838                                         fHistoMCAllGammaPt[fiCut]->Fill(particle->Pt()); // All MC Gamma
1839                                         if(particle->GetMother(0) >-1){
1840                                                 if (fMCStack->Particle(particle->GetMother(0))->GetPdgCode() ==111){
1841                                                         if (fMCStack->Particle(particle->GetMother(0))->GetMother(0) > -1){
1842                                                                 if ( fMCStack->Particle((fMCStack->Particle(particle->GetMother(0)))->GetMother(0))->GetPdgCode() == 221 ||
1843                                                                         fMCStack->Particle((fMCStack->Particle(particle->GetMother(0)))->GetMother(0))->GetPdgCode() == 223 ){
1844                                                                         if ( fMCStack->Particle(particle->GetMother(0))->GetNDaughters()==3 )
1845                                                                                 fHistoMCGammaFromNeutralMesonPt[fiCut]->Fill(particle->Pt()); // All photons from eta or omega via pi0
1846                                                                 }
1847                                                         }
1848                                                 }
1849                                         }
1850                                 }
1851                         }
1852                         if (fNeutralPionMode < 2){
1853                                 if (((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->PhotonIsSelectedMC(particle,fMCStack,kTRUE)){
1854                                         fHistoMCConvGammaPt[fiCut]->Fill(particle->Pt());
1855                                 } // Converted MC Gamma
1856                         }
1857                         if(((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->PionIsSelectedMC(i,fMCStack)){
1858                                 if( particle->GetPdgCode() == 211){
1859                                         fHistoMCAllPosPionsPt[fiCut]->Fill(particle->Pt()); // All pos pions
1860                                         if(particle->GetMother(0) >-1){
1861                                                 if (fMCStack->Particle(particle->GetMother(0))->GetPdgCode() ==221 || fMCStack->Particle(particle->GetMother(0))->GetPdgCode() ==223)
1862                                                         fHistoMCPosPionsFromNeutralMesonPt[fiCut]->Fill(particle->Pt()); // All pos from eta or omega
1863                                         }
1864                                 }
1865                                 if( particle->GetPdgCode() == -211){
1866                                         fHistoMCAllNegPionsPt[fiCut]->Fill(particle->Pt()); // All neg pions
1867                                         if(particle->GetMother(0) >-1){
1868                                                 if (fMCStack->Particle(particle->GetMother(0))->GetPdgCode() ==221 || fMCStack->Particle(particle->GetMother(0))->GetPdgCode() ==223 )
1869                                                         fHistoMCNegPionsFromNeutralMesonPt[fiCut]->Fill(particle->Pt()); // All pos from eta or omega
1870                                         }
1871                                 }
1872                         }
1873
1874
1875                         // \eta -> pi+ pi- \gamma
1876                         Int_t labelNeutPion = -1;
1877                         Int_t labelNegPion = -1;
1878                         Int_t labelPosPion = -1;
1879
1880                         if( ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelectedMCPiPlPiMiPiZero(particle,fMCStack,labelNegPion,labelPosPion,labelNeutPion,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift())){
1881                                 Float_t weighted= 1;
1882                                 if( ((AliPrimaryPionCuts*) fPionCutArray->At(fiCut))->DoWeights() ) {
1883                                         if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack,fInputEvent)){
1884                                                 if (particle->Pt()>0.005){
1885                                                         weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, fMCStack,fInputEvent);
1886                                                 }
1887                                         }
1888                                 }
1889                                 if(particle->GetPdgCode() == 221)fHistoMCEtaPiPlPiMiPiZeroPt[fiCut]->Fill(particle->Pt(), weighted);                                            // All MC Eta in respective decay channel
1890                                 if(particle->GetPdgCode() == 223)fHistoMCOmegaPiPlPiMiPiZeroPt[fiCut]->Fill(particle->Pt(), weighted);                                          // All MC Omega in respective decay channel
1891
1892                                 TParticle *neutPion    = fMCStack->Particle(labelNeutPion);
1893                                 TParticle *gamma1 = fMCStack->Particle(neutPion->GetDaughter(0));
1894                                 TParticle *gamma2 = fMCStack->Particle(neutPion->GetDaughter(1));
1895                                 if (fNeutralPionMode < 2){
1896                                         if(
1897                                                 ((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->PhotonIsSelectedMC(gamma1,fMCStack,kFALSE) &&                                    // test first daugther of pi0
1898                                                 ((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->PhotonIsSelectedMC(gamma2,fMCStack,kFALSE) &&                                    // test second daughter of pi0
1899                                                 ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->PionIsSelectedMC(labelNegPion,fMCStack) &&                                                             // test negative pion
1900                                                 ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->PionIsSelectedMC(labelPosPion,fMCStack)                                                                // test positive pion
1901                                         ) {
1902                                                         if(particle->GetPdgCode() == 221) fHistoMCEtaPiPlPiMiPiZeroInAccPt[fiCut]->Fill(particle->Pt(), weighted );             // MC Eta pi+ pi- pi0 with gamma's and e+e- in acc
1903                                                         if(particle->GetPdgCode() == 223) fHistoMCOmegaPiPlPiMiPiZeroInAccPt[fiCut]->Fill(particle->Pt(), weighted );           // MC Omega pi+ pi- pi0 with gamma's and e+e- in acc
1904                                         }
1905                                 } else if (fNeutralPionMode == 2){
1906                                         if(
1907                                                 ((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->ClusterIsSelectedMC(gamma1,fMCStack) &&                                      // test first daugther of pi0
1908                                                 ((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->ClusterIsSelectedMC(gamma2,fMCStack) &&                                      // test second daughter of pi0
1909                                                 ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->PionIsSelectedMC(labelNegPion,fMCStack) &&                                                             // test negative pion
1910                                                 ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->PionIsSelectedMC(labelPosPion,fMCStack)                                                                // test positive pion
1911                                         ) {
1912                                                         if(particle->GetPdgCode() == 221) fHistoMCEtaPiPlPiMiPiZeroInAccPt[fiCut]->Fill(particle->Pt(), weighted );             // MC Eta pi+ pi- pi0 with gamma's and e+e- in acc
1913                                                         if(particle->GetPdgCode() == 223) fHistoMCOmegaPiPlPiMiPiZeroInAccPt[fiCut]->Fill(particle->Pt(), weighted );           // MC Omega pi+ pi- pi0 with gamma's and e+e- in acc
1914                                         }
1915                                 }
1916
1917                         }
1918                 }
1919         }
1920 }
1921
1922
1923 //________________________________________________________________________
1924 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::CalculateMesonCandidates(){
1925
1926         // Conversion Gammas
1927         if( fNeutralPionCandidates->GetEntries() > 0 && fGoodVirtualParticles->GetEntries() > 0 ){
1928
1929                 for(Int_t mesonIndex=0; mesonIndex<fNeutralPionCandidates->GetEntries(); mesonIndex++){
1930                         AliAODConversionMother *neutralPion=dynamic_cast<AliAODConversionMother*>(fNeutralPionCandidates->At(mesonIndex));
1931                         if (neutralPion==NULL) continue;
1932                         for(Int_t virtualParticleIndex=0;virtualParticleIndex<fGoodVirtualParticles->GetEntries();virtualParticleIndex++){
1933
1934                                 AliAODConversionPhoton *vParticle=dynamic_cast<AliAODConversionPhoton*>(fGoodVirtualParticles->At(virtualParticleIndex));
1935                                 if (vParticle==NULL) continue;
1936                                 //Check for same Electron ID
1937
1938                                 AliAODConversionMother *mesoncand = new AliAODConversionMother(neutralPion,vParticle);
1939                                 mesoncand->SetLabels(mesonIndex,virtualParticleIndex);
1940
1941                                 if( ( ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(mesoncand,kTRUE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift())) ){
1942
1943 //                                      cout<< "Meson Accepted "<<endl;
1944
1945                                         Int_t zbin= fBGHandler[fiCut]->GetZBinIndex(fESDEvent->GetPrimaryVertex()->GetZ());
1946                                         Int_t mbin = 0;
1947                                         if( ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
1948                                                 mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fNumberOfESDTracks);
1949                                         } else {
1950                                                 if (fNeutralPionMode < 2) mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGoodConvGammas->GetEntries());
1951                                                 else mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fClusterCandidates->GetEntries());
1952                                         }
1953
1954                                         fHistoMotherInvMassPt[fiCut]->Fill(mesoncand->M(),mesoncand->Pt());
1955                                         Double_t sparesFill[4] = {mesoncand->M(),mesoncand->Pt(),(Double_t)zbin,(Double_t)mbin};
1956                                         fTHnSparseMotherInvMassPtZM[fiCut]->Fill(sparesFill,1);
1957
1958                                         if(fMCEvent){
1959                                                 ProcessTrueMesonCandidates(mesoncand,neutralPion,vParticle);
1960                                         }
1961                                 }
1962                                 delete mesoncand;
1963                                 mesoncand=0x0;
1964                         }
1965                 }
1966         }
1967 }
1968
1969 //________________________________________________________________________
1970 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::CalculateBackground(){
1971
1972         Int_t zbin= fBGHandler[fiCut]->GetZBinIndex(fESDEvent->GetPrimaryVertex()->GetZ());
1973         Int_t mbin = 0;
1974
1975
1976         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
1977                 mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fNumberOfESDTracks);
1978         } else {
1979                 if (fNeutralPionMode < 2) mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGoodConvGammas->GetEntries());
1980                 else mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fClusterCandidates->GetEntries());
1981         }
1982
1983         Int_t method = 1;
1984         AliGammaConversionAODBGHandler::GammaConversionVertex *bgEventVertex = NULL;
1985         if( ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity() ) {
1986                 for(Int_t nEventsInBG=0;nEventsInBG<fBGHandler[fiCut]->GetNBGEvents();nEventsInBG++){
1987                         AliGammaConversionMotherAODVector *previousEventMesons = fBGHandler[fiCut]->GetBGGoodMesons(zbin,mbin,nEventsInBG);
1988                         if(fMoveParticleAccordingToVertex == kTRUE && method == 1){
1989                                 bgEventVertex = fBGHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG);
1990                         }
1991
1992                         for(Int_t iCurrent=0;iCurrent<fGoodVirtualParticles->GetEntries();iCurrent++){
1993                                 AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGoodVirtualParticles->At(iCurrent));
1994
1995                                 for(UInt_t iPrevious=0;iPrevious<previousEventMesons->size();iPrevious++){
1996                                         AliAODConversionMother previousGoodMeson = (AliAODConversionMother)(*(previousEventMesons->at(iPrevious)));
1997
1998                                         if(fMoveParticleAccordingToVertex == kTRUE && method == 1 ){
1999                                                 MoveParticleAccordingToVertex(&previousGoodMeson,bgEventVertex);
2000                                         }
2001
2002                                         AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&previousGoodMeson,&currentEventGoodV0);
2003
2004
2005                                         if( ( ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE, ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift()))){
2006                                                 fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
2007                                                 Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin};
2008                                                 fTHnSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
2009                                         }
2010                                         delete backgroundCandidate;
2011                                         backgroundCandidate = 0x0;
2012                                 }
2013                         }
2014                 }
2015         } else {
2016                 for(Int_t nEventsInBG=0;nEventsInBG <fBGHandler[fiCut]->GetNBGEvents();nEventsInBG++){
2017                         AliGammaConversionMotherAODVector *previousEventMesons = fBGHandler[fiCut]->GetBGGoodMesons(zbin,mbin,nEventsInBG);
2018                         if(previousEventMesons){
2019                                 if(fMoveParticleAccordingToVertex == kTRUE && method == 1){
2020                                         bgEventVertex = fBGHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG);
2021                                 }
2022                                 for(Int_t iCurrent=0;iCurrent<fGoodVirtualParticles->GetEntries();iCurrent++){
2023
2024                                         AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGoodVirtualParticles->At(iCurrent));
2025
2026                                         for(UInt_t iPrevious=0;iPrevious<previousEventMesons->size();iPrevious++){
2027
2028                                                 AliAODConversionMother previousGoodMeson = (AliAODConversionMother)(*(previousEventMesons->at(iPrevious)));
2029
2030                                                 if(fMoveParticleAccordingToVertex == kTRUE && method ==1){
2031                                                         MoveParticleAccordingToVertex(&previousGoodMeson,bgEventVertex);
2032                                                 }
2033
2034                                                 AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&previousGoodMeson,&currentEventGoodV0);
2035
2036                                                 if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift()))){
2037                                                         fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
2038                                                         Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin};
2039                                                         fTHnSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
2040                                                 }
2041                                                 delete backgroundCandidate;
2042                                                 backgroundCandidate = 0x0;
2043                                         }
2044                                 }
2045                         }
2046                 }
2047         }
2048 }
2049
2050 //______________________________________________________________________
2051 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessTrueMesonCandidates(AliAODConversionMother *mesoncand, AliAODConversionMother *TrueNeutralPionCandidate, AliAODConversionPhoton *TrueVirtualParticleCandidate){
2052
2053         // Process True Mesons
2054         AliStack *MCStack = fMCEvent->Stack();
2055
2056         Bool_t isTrueEta = kFALSE;
2057         Bool_t isTrueOmega = kFALSE;
2058         Int_t trueMesonFlag = TrueNeutralPionCandidate->GetTrueMesonValue();
2059         Int_t pi0MCLabel= TrueNeutralPionCandidate->GetMCLabel();
2060
2061
2062         if ( !(trueMesonFlag == 1 && pi0MCLabel != -1)) return;
2063 //      cout << trueMesonFlag << "\t" << pi0MCLabel << endl;
2064
2065
2066         Int_t virtualParticleMCLabel = TrueVirtualParticleCandidate->GetMCParticleLabel(MCStack);
2067         Int_t virtualParticleMotherLabel = -1;
2068         Bool_t isPiPiDecay = kFALSE;
2069
2070 //      if (fDoMesonQA){
2071                 TParticle * negativeMC = (TParticle*)TrueVirtualParticleCandidate->GetNegativeMCDaughter(MCStack);
2072                 TParticle * positiveMC = (TParticle*)TrueVirtualParticleCandidate->GetPositiveMCDaughter(MCStack);
2073 //      }
2074
2075         if(virtualParticleMCLabel != -1){ // if virtualParticleMCLabel==-1 particles don't have same mother
2076 //              TParticle * negativeMC = (TParticle*)TrueVirtualParticleCandidate->GetNegativeMCDaughter(MCStack);
2077 //              TParticle * positiveMC = (TParticle*)TrueVirtualParticleCandidate->GetPositiveMCDaughter(MCStack);
2078 //              TParticle * virtualParticleMotherMC = (TParticle*)MCStack->Particle(virtualParticleMCLabel);
2079 //              cout << "pdg code same mother - " << virtualParticleMotherMC->GetPdgCode() << endl;
2080
2081                 if(TMath::Abs(negativeMC->GetPdgCode())==211 && TMath::Abs(positiveMC->GetPdgCode())==211){  // Pions ...
2082                         virtualParticleMotherLabel=virtualParticleMCLabel;
2083                         isPiPiDecay=kTRUE;
2084 //                      } else if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){  // Electrons ...
2085 //                              if( virtualParticleMotherMC->GetPdgCode() != 22 ){
2086 //                                      virtualParticleMotherLabel=virtualParticleMCLabel;
2087 //                                      isDalitz = kTRUE;
2088 //                              } else if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
2089 //                                      virtualParticleMotherLabel=virtualParticleMotherMC->GetFirstMother();
2090 //                                      isRealGamma = kTRUE; //no virtual gamma
2091 //                              }
2092                 }
2093         }
2094         if (IsEtaPiPlPiMiPiZeroDaughter(pi0MCLabel) || IsOmegaPiPlPiMiPiZeroDaughter(pi0MCLabel)){
2095                 Int_t pi0MotherMCLabel = ((TParticle*)MCStack->Particle(pi0MCLabel))->GetMother(0);
2096                 if(virtualParticleMCLabel != -1){
2097 //                      cout << "pi+pi- mother: "<<  virtualParticleMCLabel << endl;
2098 //                      cout << "pi0 mother: "<<  pi0MotherMCLabel << endl;
2099
2100 //                      TParticle * virtualParticleMotherMC = (TParticle*)MCStack->Particle(virtualParticleMCLabel);
2101 //                      cout << "pdg code same mother - " << virtualParticleMotherMC->GetPdgCode() << endl;
2102                 }
2103                 if( pi0MotherMCLabel == virtualParticleMotherLabel ){
2104                         if(((TParticle*)MCStack->Particle(virtualParticleMotherLabel))->GetPdgCode() == 221){
2105 //                              cout << "found eta" << endl;
2106                                 isTrueEta=kTRUE;
2107                         }
2108                         if(((TParticle*)MCStack->Particle(virtualParticleMotherLabel))->GetPdgCode() == 223){
2109 //                              cout << "found omega" << endl;
2110                                 isTrueOmega=kTRUE;
2111                         }
2112                 }
2113         }
2114
2115
2116         if( isTrueEta || isTrueOmega ){ // True Eta or Omega
2117                 if ( isPiPiDecay) { //real eta -> Pi+ Pi- Pi0
2118                         Float_t weighted= 1;
2119 //                      if( ((AliPrimaryPionCuts*) fPionCutArray->At(fiCut))->DoWeights() ) {
2120 //                              if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(gammaMotherLabel, fMCStack,fInputEvent)){
2121 //                                      if (((TParticle*)MCStack->Particle(gammaMotherLabel))->Pt()>0.005){
2122 //                                              weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),gammaMotherLabel,fMCStack,fInputEvent);
2123 //                                      }
2124 //                              }
2125 //                      }
2126                         fHistoTrueMotherPiPlPiMiPiZeroInvMassPt[fiCut]->Fill(mesoncand->M(),mesoncand->Pt(),weighted);
2127                 }
2128         }
2129
2130 }
2131
2132
2133 //________________________________________________________________________
2134 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::UpdateEventByEventData(){
2135    //see header file for documentation
2136
2137         Int_t method = 1;
2138         if( method == 1 ) {
2139                 if(fNeutralPionCandidates->GetEntries() >0 ){
2140                         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2141                                 fBGHandler[fiCut]->AddMesonEvent(fNeutralPionCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fV0Reader->GetNumberOfPrimaryTracks(),0);
2142                         } else { // means we use #V0s for multiplicity
2143                                 if (fNeutralPionMode < 2) fBGHandler[fiCut]->AddMesonEvent(fNeutralPionCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fGoodConvGammas->GetEntries(),0);
2144                                 else fBGHandler[fiCut]->AddMesonEvent(fNeutralPionCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fClusterCandidates->GetEntries(),0);
2145                         }
2146                 }
2147         } else if ( method == 2 ){
2148                 if(fGoodVirtualParticles->GetEntries() > 0 ){
2149                         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2150                                 fBGHandler[fiCut]->AddEvent(fGoodVirtualParticles,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fV0Reader->GetNumberOfPrimaryTracks(),0);
2151                         } else{ // means we use #V0s for multiplicity
2152                                 fBGHandler[fiCut]->AddEvent(fGoodVirtualParticles,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fGoodVirtualParticles->GetEntries(),0);
2153                         }
2154                 }
2155         }
2156 }
2157
2158 //________________________________________________________________________
2159 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::MoveParticleAccordingToVertex(AliAODConversionMother* particle,const AliGammaConversionAODBGHandler::GammaConversionVertex *vertex){
2160    //see header file for documentation
2161
2162    Double_t dx = vertex->fX - fInputEvent->GetPrimaryVertex()->GetX();
2163    Double_t dy = vertex->fY - fInputEvent->GetPrimaryVertex()->GetY();
2164    Double_t dz = vertex->fZ - fInputEvent->GetPrimaryVertex()->GetZ();
2165
2166    Double_t movedPlace[3] = {particle->GetProductionX() - dx,particle->GetProductionY() - dy,particle->GetProductionZ() - dz};
2167    particle->SetProductionPoint(movedPlace);
2168 }
2169
2170 //_____________________________________________________________________________________
2171 Bool_t AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::IsEtaPiPlPiMiPiZeroDaughter( Int_t label ) const {
2172 //
2173 // Returns true if the particle comes from eta -> pi+ pi- gamma
2174 //
2175         Int_t motherLabel = fMCStack->Particle( label )->GetMother(0);
2176         if( motherLabel < 0 || motherLabel >= fMCStack->GetNtrack() ) return kFALSE;
2177
2178         TParticle* mother = fMCStack->Particle( motherLabel );
2179 //      cout << "found eta? " << endl;
2180         if( mother->GetPdgCode() != 221 ) return kFALSE;
2181 //              else cout << "YES" << endl;
2182         if( IsPiPlPiMiPiZeroDecay( mother ) ) return kTRUE;
2183         return kFALSE;
2184 }
2185
2186 //_____________________________________________________________________________________
2187 Bool_t AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::IsOmegaPiPlPiMiPiZeroDaughter( Int_t label ) const {
2188 //
2189 // Returns true if the particle comes from eta -> pi+ pi- gamma
2190 //
2191         Int_t motherLabel = fMCStack->Particle( label )->GetMother(0);
2192         if( motherLabel < 0 || motherLabel >= fMCStack->GetNtrack() ) return kFALSE;
2193
2194         TParticle* mother = fMCStack->Particle( motherLabel );
2195 //      cout << "found omega? " << endl;
2196         if( mother->GetPdgCode() != 223 ) return kFALSE;
2197 //              else cout << "YES" << endl;
2198         if( IsPiPlPiMiPiZeroDecay( mother ) ) return kTRUE;
2199         return kFALSE;
2200 }
2201
2202
2203 //_____________________________________________________________________________
2204 Bool_t AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::IsPiPlPiMiPiZeroDecay(TParticle *fMCMother) const
2205 {
2206 //      cout << fMCMother->GetNDaughters() << endl;
2207         if( fMCMother->GetNDaughters() != 3 ) return kFALSE;
2208 //      cout << fMCMother->GetPdgCode() << endl;
2209         if( !(fMCMother->GetPdgCode() == 221 || fMCMother->GetPdgCode() == 223)  ) return kFALSE;
2210 //      cout << "made it til here" << endl;
2211
2212         TParticle *posPion = 0x0;
2213         TParticle *negPion = 0x0;
2214         TParticle *neutPion    = 0x0;
2215
2216         for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
2217                 TParticle* temp = (TParticle*)fMCStack->Particle( index );
2218
2219                 switch( temp->GetPdgCode() ) {
2220                 case 211:
2221                         posPion =  temp;
2222                         break;
2223                 case -211:
2224                         negPion =  temp;
2225                         break;
2226                 case 111:
2227                         neutPion = temp;
2228                         break;
2229                 }
2230         }
2231         if( posPion && negPion && neutPion) return kTRUE;
2232
2233         return kFALSE;
2234 }
2235
2236 //_____________________________________________________________________________________
2237 Bool_t AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::GammaIsNeutralMesonPiPlPiMiPiZeroDaughter( Int_t label ) const {
2238 //
2239 // Returns true if the particle comes from eta -> pi+ pi- gamma
2240 //
2241         Int_t motherLabel = fMCStack->Particle( label )->GetMother(0);
2242         if( motherLabel < 0 || motherLabel >= fMCStack->GetNtrack() ) return kFALSE;
2243
2244         TParticle* mother = fMCStack->Particle( motherLabel );
2245 //      cout << "found omega? " << endl;
2246         if( mother->GetPdgCode() != 111 ) return kFALSE;
2247 //              else cout << "YES" << endl;
2248         Int_t grandMotherLabel = mother->GetMother(0);
2249         if( grandMotherLabel < 0 || grandMotherLabel >= fMCStack->GetNtrack() ) return kFALSE;
2250         TParticle* grandmother = fMCStack->Particle( grandMotherLabel );
2251
2252         if( IsPiPlPiMiPiZeroDecay( grandmother ) ) return kTRUE;
2253         return kFALSE;
2254 }