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