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
  20 #include "TParticle.h"
  21 #include "TPDGCode.h"
  22 #include "TMCProcess.h"
  23 #include "TDatabasePDG.h"
  24 #include "TList.h"
  25 #include "TChain.h"
  26 #include "TDirectory.h"
  27 #include "TTree.h"
  28 #include "TH1.h"
  29 #include "TH1F.h"
  30 #include "THnSparse.h"
  31 #include "TH2F.h"
  32 #include "AliStack.h"
  33 #include "AliAnalysisManager.h"
  34 #include "AliESDInputHandler.h"
  35 #include "AliESDtrack.h"
  36 #include "AliMCEvent.h"
  37 #include "AliStack.h"
  38 #include "AliMCEventHandler.h"
  39 #include "AliPID.h"
  40 #include "AliLog.h"
  41 #include "AliESDtrackCuts.h"
  42 #include "AliESDpidCuts.h"
  43 #include "AliMCEvent.h"
  44 #include "AliESDv0.h"
  45 #include "AliESDEvent.h"
  46 #include "AliESDpid.h"
  47 #include "AliKFParticle.h"
  48 #include "AliMCEventHandler.h"
  49 #include "AliKFVertex.h"
  50 #include "AliTriggerAnalysis.h"
  51 #include "AliCentrality.h"
  52 #include "AliMultiplicity.h"
  53 #include "AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero.h"
  54
  55
  56 ClassImp( AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero )
  57
  58 //-----------------------------------------------------------------------------------------------
  59 AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero():
  60         fV0Reader(NULL),
  61         fPionSelector(NULL),
  62         fBGHandler(NULL),
  63         fESDEvent(NULL),
  64         fMCEvent(NULL),
  65         fMCStack(NULL),
  66         fCutFolder(NULL),
  67         fESDList(NULL),
  68         fBackList(NULL),
  69         fMotherList(NULL),
  70         fTrueList(NULL),
  71         fMCList(NULL),
  72         fOutputContainer(0),
  73         fReaderGammas(NULL),
  74         fSelectorNegPionIndex(0),
  75         fSelectorPosPionIndex(0),
  76         fGoodGammas(NULL),
  77         fNeutralPionCandidates(NULL),
  78         fGoodVirtualParticles(NULL),
  79         fEventCutArray(NULL),
  80         fGammaCutArray(NULL),
  81         fPionCutArray(NULL),
  82         fNeutralPionMesonCutArray(NULL),
  83         fMesonCutArray(NULL),
  84         fEventCuts(NULL),
  85         fConversionCuts(NULL),
  86         fHistoConvGammaPt(NULL),
  87         fHistoConvGammaEta(NULL),
  88         fHistoNegPionPt(NULL),
  89         fHistoPosPionPt(NULL),
  90         fHistoNegPionPhi(NULL),
  91         fHistoPosPionPhi(NULL),
  92         fHistoNegPionEta(NULL),
  93         fHistoPosPionEta(NULL),
  94         fHistoNegPionClsTPC(NULL),
  95         fHistoPosPionClsTPC(NULL),
  96         fHistoPionDCAxy(NULL),
  97         fHistoPionDCAz(NULL),
  98         fHistoPionTPCdEdxNSigma(NULL),
  99         fHistoPionTPCdEdx(NULL),
 100         fHistoPionPionInvMassPt(NULL),
 101         fHistoGammaGammaInvMassPt(NULL),
 102         fHistoMotherInvMassPt(NULL),
 103         fTHnSparseMotherInvMassPtZM(NULL),
 104         fHistoMotherBackInvMassPt(NULL),
 105         fTHnSparseMotherBackInvMassPtZM(NULL),
 106         fHistoMCAllGammaPt(NULL),
 107         fHistoMCConvGammaPt(NULL),
 108         fHistoMCAllPosPionsPt(NULL),
 109         fHistoMCAllNegPionsPt(NULL),
 110         fHistoMCGammaFromNeutralMesonPt(NULL),
 111         fHistoMCPosPionsFromNeutralMesonPt(NULL),
 112         fHistoMCNegPionsFromNeutralMesonPt(NULL),
 113         fHistoMCEtaPiPlPiMiPiZeroPt(NULL),
 114         fHistoMCEtaPiPlPiMiPiZeroInAccPt(NULL),
 115         fHistoMCOmegaPiPlPiMiPiZeroPt(NULL),
 116         fHistoMCOmegaPiPlPiMiPiZeroInAccPt(NULL),
 117         fHistoTrueMotherPiPlPiMiPiZeroInvMassPt(NULL),
 118         fHistoTrueMotherGammaGammaInvMassPt(NULL),
 119         fHistoTrueConvGammaPt(NULL),
 120         fHistoTrueConvGammaFromNeutralMesonPt(NULL),
 121         fHistoTruePosPionPt(NULL),
 122         fHistoTruePosPionFromNeutralMesonPt(NULL),
 123         fHistoTrueNegPionPt(NULL),
 124         fHistoTrueNegPionFromNeutralMesonPt(NULL),
 125         fHistoTruePionPionInvMassPt(NULL),
 126         fHistoTruePionPionFromNeutralMesonInvMassPt(NULL),
 127         fHistoNEvents(NULL),
 128         fHistoNGoodESDTracks(NULL),
 129         fProfileEtaShift(NULL),
 130         fRandom(0),
 131         fnCuts(0),
 132         fiCut(0),
 133         fNumberOfESDTracks(0),
 134         fMoveParticleAccordingToVertex(kFALSE),
 135         fIsHeavyIon(0),
 136         fDoMesonAnalysis(kTRUE),
 137         fDoMesonQA(kFALSE),
 138         fIsFromMBHeader(kTRUE),
 139         fIsMC(kFALSE)
 140 {
 141
 142 }
 143
 144 //-----------------------------------------------------------------------------------------------
 145 AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero( const char* name ):
 146         AliAnalysisTaskSE(name),
 147         fV0Reader(NULL),
 148         fPionSelector(NULL),
 149         fBGHandler(NULL),
 150         fESDEvent(NULL),
 151         fMCEvent(NULL),
 152         fMCStack(NULL),
 153         fCutFolder(NULL),
 154         fESDList(NULL),
 155         fBackList(NULL),
 156         fMotherList(NULL),
 157         fTrueList(NULL),
 158         fMCList(NULL),
 159         fOutputContainer(0),
 160         fReaderGammas(NULL),
 161         fSelectorNegPionIndex(0),
 162         fSelectorPosPionIndex(0),
 163         fGoodGammas(NULL),
 164         fNeutralPionCandidates(NULL),
 165         fGoodVirtualParticles(NULL),
 166         fEventCutArray(NULL),
 167         fGammaCutArray(NULL),
 168         fPionCutArray(NULL),
 169         fNeutralPionMesonCutArray(NULL),
 170         fMesonCutArray(NULL),
 171         fEventCuts(NULL),
 172         fConversionCuts(NULL),
 173         fHistoConvGammaPt(NULL),
 174         fHistoConvGammaEta(NULL),
 175         fHistoNegPionPt(NULL),
 176         fHistoPosPionPt(NULL),
 177         fHistoNegPionPhi(NULL),
 178         fHistoPosPionPhi(NULL),
 179         fHistoNegPionEta(NULL),
 180         fHistoPosPionEta(NULL),
 181         fHistoNegPionClsTPC(NULL),
 182         fHistoPosPionClsTPC(NULL),
 183         fHistoPionDCAxy(NULL),
 184         fHistoPionDCAz(NULL),
 185         fHistoPionTPCdEdxNSigma(NULL),
 186         fHistoPionTPCdEdx(NULL),
 187         fHistoPionPionInvMassPt(NULL),
 188         fHistoGammaGammaInvMassPt(NULL),
 189         fHistoMotherInvMassPt(NULL),
 190         fTHnSparseMotherInvMassPtZM(NULL),
 191         fHistoMotherBackInvMassPt(NULL),
 192         fTHnSparseMotherBackInvMassPtZM(NULL),
 193         fHistoMCAllGammaPt(NULL),
 194         fHistoMCConvGammaPt(NULL),
 195         fHistoMCAllPosPionsPt(NULL),
 196         fHistoMCAllNegPionsPt(NULL),
 197         fHistoMCGammaFromNeutralMesonPt(NULL),
 198         fHistoMCPosPionsFromNeutralMesonPt(NULL),
 199         fHistoMCNegPionsFromNeutralMesonPt(NULL),
 200         fHistoMCEtaPiPlPiMiPiZeroPt(NULL),
 201         fHistoMCEtaPiPlPiMiPiZeroInAccPt(NULL),
 202         fHistoMCOmegaPiPlPiMiPiZeroPt(NULL),
 203         fHistoMCOmegaPiPlPiMiPiZeroInAccPt(NULL),
 204         fHistoTrueMotherPiPlPiMiPiZeroInvMassPt(NULL),
 205         fHistoTrueMotherGammaGammaInvMassPt(NULL),
 206         fHistoTrueConvGammaPt(NULL),
 207         fHistoTrueConvGammaFromNeutralMesonPt(NULL),
 208         fHistoTruePosPionPt(NULL),
 209         fHistoTruePosPionFromNeutralMesonPt(NULL),
 210         fHistoTrueNegPionPt(NULL),
 211         fHistoTrueNegPionFromNeutralMesonPt(NULL),
 212         fHistoTruePionPionInvMassPt(NULL),
 213         fHistoTruePionPionFromNeutralMesonInvMassPt(NULL),
 214         fHistoNEvents(NULL),
 215         fHistoNGoodESDTracks(NULL),
 216         fProfileEtaShift(NULL),
 217         fRandom(0),
 218         fnCuts(0),
 219         fiCut(0),
 220         fNumberOfESDTracks(0),
 221         fMoveParticleAccordingToVertex(kFALSE),
 222         fIsHeavyIon(0),
 223         fDoMesonAnalysis(kTRUE),
 224         fDoMesonQA(kFALSE),
 225         fIsFromMBHeader(kTRUE),
 226         fIsMC(kFALSE)
 227 {
 228    DefineOutput(1, TList::Class());
 229 }
 230
 231 //-----------------------------------------------------------------------------------------------
 232 AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::~AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero()
 233 {
 234         //
 235         // virtual destructor
 236         //
 237         cout<<"Destructor"<<endl;
 238         if(fGoodGammas){
 239                 delete fGoodGammas;
 240                 fGoodGammas = 0x0;
 241         }
 242         if(fNeutralPionCandidates){
 243                 delete fNeutralPionCandidates;
 244                 fNeutralPionCandidates = 0x0;
 245         }
 246
 247         if(fGoodVirtualParticles){
 248                 delete fGoodVirtualParticles;
 249                 fGoodVirtualParticles = 0x0;
 250         }
 251         if(fBGHandler){
 252                 delete[] fBGHandler;
 253                 fBGHandler = 0x0;
 254         }
 255 }
 256 //___________________________________________________________
 257 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::InitBack(){
 258
 259         const Int_t nDim = 4;
 260         Int_t nBins[nDim] = {500,250,7,4};
 261         Double_t xMin[nDim] = {0.4,0, 0,0};
 262         Double_t xMax[nDim] = {0.9,25,7,4};
 263
 264         fTHnSparseMotherInvMassPtZM = new THnSparseF*[fnCuts];
 265         fTHnSparseMotherBackInvMassPtZM = new THnSparseF*[fnCuts];
 266
 267         fBGHandler = new AliGammaConversionAODBGHandler*[fnCuts];
 268         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 269
 270                 TString cutstringEvent          = ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber();
 271                 TString cutstringPion           = ((AliPrimaryPionCuts*)fPionCutArray->At(iCut))->GetCutNumber();
 272                 TString cutstringGamma          = ((AliConversionPhotonCuts*)fGammaCutArray->At(iCut))->GetCutNumber();
 273                 TString cutstringNeutralPion= ((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(iCut))->GetCutNumber();
 274                 TString cutstringMeson          = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
 275
 276                 TString fullCutString = Form("%s_%s_%s_%s_%s",cutstringEvent.Data(),cutstringGamma.Data(),cutstringNeutralPion.Data(), cutstringPion.Data(),cutstringMeson.Data());
 277                 TString nameBackList = Form("%s Back histograms",fullCutString.Data());
 278                 TString nameMotherList = Form("%s Mother histograms",fullCutString.Data());
 279
 280                 Int_t collisionSystem = atoi((TString)(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber())(0,1));
 281                 Int_t centMin = atoi((TString)(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber())(1,1));
 282                 Int_t centMax = atoi((TString)(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber())(2,1));
 283
 284                 if(collisionSystem == 1 || collisionSystem == 2 ||
 285                         collisionSystem == 5 || collisionSystem == 8 ||
 286                         collisionSystem == 9){
 287                         centMin = centMin*10;
 288                         centMax = centMax*10;
 289                 }
 290                 else if(collisionSystem == 3 || collisionSystem == 6){
 291                         centMin = centMin*5;
 292                         centMax = centMax*5;
 293                 }
 294                 else if(collisionSystem == 4 || collisionSystem == 7){
 295                         centMin = ((centMin*5)+45);
 296                         centMax = ((centMax*5)+45);
 297                 }
 298
 299
 300                 fBackList[iCut] = new TList();
 301
 302                 fBackList[iCut]->SetName(nameBackList.Data());
 303                 fBackList[iCut]->SetOwner(kTRUE);
 304                 fCutFolder[iCut]->Add(fBackList[iCut]);
 305
 306                 fTHnSparseMotherBackInvMassPtZM[iCut] = new THnSparseF("Back_Back_InvMass_Pt_z_m","Back_Back_InvMass_Pt_z_m",nDim,nBins,xMin,xMax);
 307                 fBackList[iCut]->Add(fTHnSparseMotherBackInvMassPtZM[iCut]);
 308
 309                 fMotherList[iCut] = new TList();
 310                 fMotherList[iCut]->SetName(nameMotherList.Data());
 311                 fMotherList[iCut]->SetOwner(kTRUE);
 312                 fCutFolder[iCut]->Add(fMotherList[iCut]);
 313
 314                 fTHnSparseMotherInvMassPtZM[iCut] = new THnSparseF("Back_Mother_InvMass_Pt_z_m","Back_Mother_InvMass_Pt_z_m",nDim,nBins,xMin,xMax);
 315                 fMotherList[iCut]->Add(fTHnSparseMotherInvMassPtZM[iCut]);
 316
 317
 318                 fBGHandler[iCut] = new AliGammaConversionAODBGHandler(  collisionSystem,centMin,centMax,
 319                                                                                                                                 ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->GetNumberOfBGEvents(),
 320                                                                                                                                 ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity());
 321         }
 322 }
 323
 324 //______________________________________________________________________
 325 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::UserCreateOutputObjects()
 326 {
 327         //
 328         // Create ouput objects
 329         //
 330
 331         // Create the output container
 332         if(fOutputContainer != NULL){
 333                 delete fOutputContainer;
 334                 fOutputContainer = NULL;
 335         }
 336         if(fOutputContainer == NULL){
 337                 fOutputContainer = new TList();
 338                 fOutputContainer->SetOwner(kTRUE);
 339         }
 340
 341         fGoodGammas = new TList();
 342         fNeutralPionCandidates = new TList();
 343         fGoodVirtualParticles = new TList();
 344
 345         fCutFolder                              = new TList*[fnCuts];
 346         fESDList                                = new TList*[fnCuts];
 347         fBackList                               = new TList*[fnCuts];
 348         fMotherList                     = new TList*[fnCuts];
 349         fHistoNEvents                   = new TH1I*[fnCuts];
 350         fHistoNGoodESDTracks    = new TH1I*[fnCuts];
 351         fProfileEtaShift                = new TProfile*[fnCuts];
 352         fHistoConvGammaPt               = new TH1F*[fnCuts];
 353         fHistoConvGammaEta              = new TH1F*[fnCuts];
 354         fHistoNegPionPt                 = new TH1F*[fnCuts];
 355         fHistoPosPionPt                 = new TH1F*[fnCuts];
 356         fHistoNegPionPhi                = new TH1F*[fnCuts];
 357         fHistoPosPionPhi                = new TH1F*[fnCuts];
 358
 359         if( fDoMesonQA ) {
 360                 fHistoNegPionEta                = new TH1F*[fnCuts];
 361                 fHistoPosPionEta                = new TH1F*[fnCuts];
 362                 fHistoNegPionClsTPC             = new TH2F*[fnCuts];
 363                 fHistoPosPionClsTPC             = new TH2F*[fnCuts];
 364                 fHistoPionDCAxy                 = new TH2F*[fnCuts];
 365                 fHistoPionDCAz                  = new TH2F*[fnCuts];
 366                 fHistoPionTPCdEdxNSigma = new TH2F*[fnCuts];
 367                 fHistoPionTPCdEdx               = new TH2F*[fnCuts];
 368                 fHistoPionPionInvMassPt = new TH2F*[fnCuts];
 369
 370         }
 371         fHistoGammaGammaInvMassPt       = new TH2F*[fnCuts];
 372         fHistoMotherInvMassPt           = new TH2F*[fnCuts];
 373         fHistoMotherBackInvMassPt       = new TH2F*[fnCuts];
 374
 375         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 376                 TString cutstringEvent          = ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber();
 377                 TString cutstringPion           = ((AliPrimaryPionCuts*)fPionCutArray->At(iCut))->GetCutNumber();
 378                 TString cutstringGamma          = ((AliConversionPhotonCuts*)fGammaCutArray->At(iCut))->GetCutNumber();
 379                 TString cutstringNeutralPion= ((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(iCut))->GetCutNumber();
 380                 TString cutstringMeson          = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
 381
 382                 TString fullCutString = Form("%s_%s_%s_%s_%s",cutstringEvent.Data(),cutstringGamma.Data(),cutstringNeutralPion.Data(), cutstringPion.Data(),cutstringMeson.Data());
 383                 TString nameCutFolder = Form("Cut Number %s", fullCutString.Data());
 384                 TString nameESDList = Form("%s ESD histograms", fullCutString.Data());
 385
 386
 387                 fCutFolder[iCut] = new TList();
 388                 fCutFolder[iCut]->SetName(nameCutFolder.Data());
 389                 fCutFolder[iCut]->SetOwner(kTRUE);
 390                 fOutputContainer->Add(fCutFolder[iCut]);
 391
 392                 fESDList[iCut] = new TList();
 393                 fESDList[iCut]->SetName(nameESDList.Data());
 394                 fESDList[iCut]->SetOwner(kTRUE);
 395
 396                 fHistoNEvents[iCut] = new TH1I("NEvents","NEvents",9,-0.5,8.5);
 397                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(1,"Accepted");
 398                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(2,"Centrality");
 399                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(3,"Missing MC");
 400                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(4,"Trigger");
 401                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(5,"Vertex Z");
 402                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(6,"Cont. Vertex");
 403                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(7,"Pile-Up");
 404                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(8,"no SDD");
 405                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(9,"no V0AND");
 406                 fESDList[iCut]->Add(fHistoNEvents[iCut]);
 407
 408                 if(fIsHeavyIon>0) fHistoNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",3000,0,3000);
 409                         else fHistoNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",200,0,200);
 410                 fESDList[iCut]->Add(fHistoNGoodESDTracks[iCut]);
 411
 412                 fProfileEtaShift[iCut] = new TProfile("Eta Shift","Eta Shift",1, -0.5,0.5);
 413                 fESDList[iCut]->Add(fProfileEtaShift[iCut]);
 414                 fHistoConvGammaPt[iCut] = new TH1F("ESD_ConvGamma_Pt","ESD_ConvGamma_Pt",250,0,25);
 415                 fESDList[iCut]->Add(fHistoConvGammaPt[iCut]);
 416                 fHistoConvGammaEta[iCut] = new TH1F("ESD_ConvGamma_Eta","ESD_ConvGamma_Eta",600,-1.5,1.5);
 417                 fESDList[iCut]->Add(fHistoConvGammaEta[iCut]);
 418                 fHistoNegPionPt[iCut] = new TH1F("ESD_PrimaryNegPions_Pt","ESD_PrimaryNegPions_Pt",1000,0,25);
 419                 fESDList[iCut]->Add(fHistoNegPionPt[iCut]);
 420                 fHistoPosPionPt[iCut] = new TH1F("ESD_PrimaryPosPions_Pt","ESD_PrimaryPosPions_Pt",1000,0,25);
 421                 fESDList[iCut]->Add(fHistoPosPionPt[iCut]);
 422                 fHistoNegPionPhi[iCut] = new TH1F("ESD_PrimaryNegPions_Phi","ESD_PrimaryNegPions_Phi",360,0,2*TMath::Pi());
 423                 fESDList[iCut]->Add(fHistoNegPionPhi[iCut]);
 424                 fHistoPosPionPhi[iCut] = new TH1F("ESD_PrimaryPosPions_Phi","ESD_PrimaryPosPions_Phi",360,0,2*TMath::Pi());
 425                 fESDList[iCut]->Add(fHistoPosPionPhi[iCut]);
 426
 427                 if ( fDoMesonQA ) {
 428                         fHistoNegPionEta[iCut] = new TH1F("ESD_PrimaryNegPions_Eta","ESD_PrimaryNegPions_Eta",600,-1.5,1.5);
 429                         fESDList[iCut]->Add(fHistoNegPionEta[iCut]);
 430                         fHistoPosPionEta[iCut] = new TH1F("ESD_PrimaryPosPions_Eta","ESD_PrimaryPosPions_Eta",600,-1.5,1.5);
 431                         fESDList[iCut]->Add(fHistoPosPionEta[iCut]);
 432                         fHistoNegPionClsTPC[iCut]  = new TH2F("ESD_PrimaryNegPions_ClsTPC","ESD_PrimaryNegPions_ClsTPC",100,0,1,400,0.,10.);
 433                         fESDList[iCut]->Add(fHistoNegPionClsTPC[iCut]);
 434                         fHistoPosPionClsTPC[iCut]  = new TH2F("ESD_PrimaryPosPions_ClsTPC","ESD_PrimaryPosPions_ClsTPC",100,0,1,400,0.,10.);
 435                         fESDList[iCut]->Add(fHistoPosPionClsTPC[iCut]);
 436                         fHistoPionDCAxy[iCut] = new TH2F("ESD_PrimaryPions_DCAxy","ESD_PrimaryPions_DCAxy",800,-4.0,4.0,400,0.,10.);
 437                         fESDList[iCut]->Add(fHistoPionDCAxy[iCut]);
 438                         fHistoPionDCAz[iCut]  = new TH2F("ESD_PrimaryPions_DCAz","ESD_PrimaryPions_DCAz",800,-4.0,4.0,400,0.,10.);
 439                         fESDList[iCut]->Add(fHistoPionDCAz[iCut]);
 440                         fHistoPionTPCdEdxNSigma[iCut] = new TH2F("ESD_PrimaryPions_TPCdEdx","ESD_PrimaryPions_TPCdEdx",150,0.05,20,400,-10,10);
 441                         fESDList[iCut]->Add(fHistoPionTPCdEdxNSigma[iCut]);
 442                         fHistoPionTPCdEdx[iCut] =new TH2F("ESD_PrimaryPions_TPCdEdxSignal","ESD_PrimaryPions_TPCdEdxSignal" ,150,0.05,20.0,800,0.0,200);
 443                         fESDList[iCut]->Add(fHistoPionTPCdEdx[iCut]);
 444                         fHistoPionPionInvMassPt[iCut] = new TH2F("ESD_PiPlusPiNeg_InvMassPt","ESD_PiPlusPiNeg_InvMassPt",2000,0.,2.,200,0.,20.);
 445                         fESDList[iCut]->Add(fHistoPionPionInvMassPt[iCut]);
 446                 }
 447                 fHistoGammaGammaInvMassPt[iCut] = new TH2F("ESD_GammaGamma_InvMass_Pt","ESD_GammaGamma_InvMass_Pt",450,0.,0.45,250,0,25);
 448                 fESDList[iCut]->Add(fHistoGammaGammaInvMassPt[iCut]);
 449                 fHistoMotherInvMassPt[iCut] = new TH2F("ESD_Mother_InvMass_Pt","ESD_Mother_InvMass_Pt",500,0.4,0.9,250,0,25);
 450                 fESDList[iCut]->Add(fHistoMotherInvMassPt[iCut]);
 451                 fHistoMotherBackInvMassPt[iCut] = new TH2F("ESD_Background_InvMass_Pt","ESD_Background_InvMass_Pt",500,0.4,0.9,250,0,25);
 452                 fESDList[iCut]->Add(fHistoMotherBackInvMassPt[iCut]);
 453
 454                 if ( fDoMesonQA ) {
 455                         TAxis *AxisAfter = fHistoPionTPCdEdxNSigma[iCut]->GetXaxis();
 456                         Int_t bins = AxisAfter->GetNbins();
 457                         Double_t from = AxisAfter->GetXmin();
 458                         Double_t to = AxisAfter->GetXmax();
 459                         Double_t *newBins = new Double_t[bins+1];
 460                         newBins[0] = from;
 461                         Double_t factor = TMath::Power(to/from, 1./bins);
 462                         for(Int_t i=1; i<=bins; ++i) newBins[i] = factor * newBins[i-1];
 463
 464                         AxisAfter->Set(bins, newBins);
 465                         AxisAfter = fHistoPionTPCdEdx[iCut]->GetXaxis();
 466                         AxisAfter->Set(bins, newBins);
 467
 468                         delete [] newBins;
 469                 }
 470
 471                 fCutFolder[iCut]->Add(fESDList[iCut]);
 472
 473         }
 474
 475         if( fIsMC ){
 476                 // MC Histogramms
 477                 fMCList = new TList*[fnCuts];
 478                 // True Histogramms
 479                 fTrueList = new TList*[fnCuts];
 480                 fHistoTrueConvGammaPt = new TH1F*[fnCuts];
 481                 fHistoTrueConvGammaFromNeutralMesonPt = new TH1F*[fnCuts];
 482                 fHistoTruePosPionPt  = new TH1F*[fnCuts];
 483                 fHistoTrueNegPionPt  = new TH1F*[fnCuts];
 484                 fHistoTruePosPionFromNeutralMesonPt  = new TH1F*[fnCuts];
 485                 fHistoTrueNegPionFromNeutralMesonPt  = new TH1F*[fnCuts];
 486
 487
 488                 fHistoMCAllGammaPt  = new TH1F*[fnCuts];
 489                 fHistoMCConvGammaPt = new TH1F*[fnCuts];
 490                 fHistoMCAllPosPionsPt = new TH1F*[fnCuts];
 491                 fHistoMCAllNegPionsPt = new TH1F*[fnCuts];
 492                 fHistoMCGammaFromNeutralMesonPt  = new TH1F*[fnCuts];
 493                 fHistoMCPosPionsFromNeutralMesonPt = new TH1F*[fnCuts];
 494                 fHistoMCNegPionsFromNeutralMesonPt = new TH1F*[fnCuts];
 495
 496 //              hMCPi0DalitzGammaPt    = new TH1F*[fnCuts];
 497 //              hMCPi0DalitzElectronPt = new TH1F*[fnCuts];
 498 //              hMCPi0DalitzPositronPt = new TH1F*[fnCuts];
 499
 500                 fHistoMCEtaPiPlPiMiPiZeroPt = new TH1F*[fnCuts];
 501                 fHistoMCEtaPiPlPiMiPiZeroInAccPt = new TH1F*[fnCuts];
 502                 fHistoMCOmegaPiPlPiMiPiZeroPt = new TH1F*[fnCuts];
 503                 fHistoMCOmegaPiPlPiMiPiZeroInAccPt = new TH1F*[fnCuts];
 504
 505                 fHistoTrueMotherPiPlPiMiPiZeroInvMassPt = new TH2F*[fnCuts];
 506                 fHistoTrueMotherGammaGammaInvMassPt = new TH2F*[fnCuts];
 507
 508 //              if (fDoMesonQA){
 509                         fHistoTruePionPionInvMassPt =                   new TH2F*[fnCuts];
 510                         fHistoTruePionPionFromNeutralMesonInvMassPt =   new TH2F*[fnCuts];
 511 //              }
 512
 513                 for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 514                         TString cutstringEvent          = ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber();
 515                         TString cutstringPion           = ((AliPrimaryPionCuts*)fPionCutArray->At(iCut))->GetCutNumber();
 516                         TString cutstringGamma          = ((AliConversionPhotonCuts*)fGammaCutArray->At(iCut))->GetCutNumber();
 517                         TString cutstringNeutralPion= ((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(iCut))->GetCutNumber();
 518                         TString cutstringMeson          = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
 519
 520                         TString fullCutString = Form("%s_%s_%s_%s_%s",cutstringEvent.Data(),cutstringGamma.Data(),cutstringNeutralPion.Data(), cutstringPion.Data(),cutstringMeson.Data());
 521                         TString nameMCList = Form("%s MC histograms", fullCutString.Data());
 522                         TString nameTrueRecList = Form("%s True histograms", fullCutString.Data());
 523
 524                         fMCList[iCut] = new TList();
 525                         fMCList[iCut]->SetName(nameMCList.Data());
 526                         fMCList[iCut]->SetOwner(kTRUE);
 527                         fCutFolder[iCut]->Add(fMCList[iCut]);
 528
 529                         fHistoMCAllGammaPt[iCut] = new TH1F("MC_AllGamma_Pt","MC_AllGamma_Pt",250,0,25);
 530                         fMCList[iCut]->Add(fHistoMCAllGammaPt[iCut]);
 531                         fHistoMCConvGammaPt[iCut] = new TH1F("MC_ConvGamma_Pt","MC_ConvGamma_Pt",250,0,25);
 532                         fMCList[iCut]->Add(fHistoMCConvGammaPt[iCut]);
 533                         fHistoMCAllPosPionsPt[iCut] = new TH1F("MC_AllPosPions_Pt","MC_AllPosPions_Pt",1000,0,25);
 534                         fMCList[iCut]->Add(fHistoMCAllPosPionsPt[iCut]);
 535                         fHistoMCAllNegPionsPt[iCut] = new TH1F("MC_AllNegPions_Pt","MC_AllNegPions_Pt",1000,0,25);
 536                         fMCList[iCut]->Add(fHistoMCAllNegPionsPt[iCut]);
 537                         fHistoMCGammaFromNeutralMesonPt[iCut] = new TH1F("MC_GammaFromNeutralMeson_Pt","MC_GammaFromNeutralMeson_Pt",250,0,25);
 538                         fMCList[iCut]->Add(fHistoMCGammaFromNeutralMesonPt[iCut]);
 539                         fHistoMCPosPionsFromNeutralMesonPt[iCut] = new TH1F("MC_PosPionsFromNeutralMeson_Pt","MC_PosPionsFromNeutralMeson_Pt",1000,0,25);
 540                         fMCList[iCut]->Add(fHistoMCPosPionsFromNeutralMesonPt[iCut]);
 541                         fHistoMCNegPionsFromNeutralMesonPt[iCut] = new TH1F("MC_NegPionsFromNeutralMeson_Pt","MC_NegPionsFromNeutralMeson_Pt",1000,0,25);
 542                         fMCList[iCut]->Add(fHistoMCNegPionsFromNeutralMesonPt[iCut]);
 543
 544                         fHistoMCEtaPiPlPiMiPiZeroPt[iCut] = new TH1F("MC_Eta_Pt","MC_Eta_Pt",250,0,25);
 545                         fHistoMCEtaPiPlPiMiPiZeroPt[iCut]->Sumw2();
 546                         fMCList[iCut]->Add(fHistoMCEtaPiPlPiMiPiZeroPt[iCut]);
 547
 548                         fHistoMCEtaPiPlPiMiPiZeroInAccPt[iCut] = new TH1F("MC_EtaInAcc_Pt","MC_EtaInAcc_Pt",250,0,25);
 549                         fHistoMCEtaPiPlPiMiPiZeroInAccPt[iCut]->Sumw2();
 550                         fMCList[iCut]->Add(fHistoMCEtaPiPlPiMiPiZeroInAccPt[iCut]);
 551
 552                         fHistoMCOmegaPiPlPiMiPiZeroPt[iCut] = new TH1F("MC_Omega_Pt","MC_Omega_Pt",250,0,25);
 553                         fHistoMCOmegaPiPlPiMiPiZeroPt[iCut]->Sumw2();
 554                         fMCList[iCut]->Add(fHistoMCOmegaPiPlPiMiPiZeroPt[iCut]);
 555
 556                         fHistoMCOmegaPiPlPiMiPiZeroInAccPt[iCut] = new TH1F("MC_OmegaInAcc_Pt","MC_OmegaInAcc_Pt",250,0,25);
 557                         fHistoMCOmegaPiPlPiMiPiZeroInAccPt[iCut]->Sumw2();
 558                         fMCList[iCut]->Add(fHistoMCOmegaPiPlPiMiPiZeroInAccPt[iCut]);
 559
 560                         fTrueList[iCut] = new TList();
 561                         fTrueList[iCut]->SetName(nameTrueRecList.Data());
 562                         fTrueList[iCut]->SetOwner(kTRUE);
 563                         fCutFolder[iCut]->Add(fTrueList[iCut]);
 564
 565                         fHistoTrueConvGammaPt[iCut] = new TH1F("ESD_TrueConvGamma_Pt","ESD_TrueConvGamma_Pt",250,0,25);
 566                         fTrueList[iCut]->Add(fHistoTrueConvGammaPt[iCut]);
 567                         fHistoTrueConvGammaFromNeutralMesonPt[iCut] = new TH1F("ESD_TrueConvGammaFromNeutralMeson_Pt","ESD_TrueConvGammaFromNeutralMeson_Pt",250,0,25);
 568                         fTrueList[iCut]->Add(fHistoTrueConvGammaFromNeutralMesonPt[iCut]);
 569
 570                         fHistoTruePosPionPt[iCut] = new TH1F("ESD_TruePosPion_Pt","ESD_TruePosPion_Pt",1000,0,25);
 571                         fTrueList[iCut]->Add(fHistoTruePosPionPt[iCut]);
 572                         fHistoTrueNegPionPt[iCut] = new TH1F("ESD_TrueNegPion_Pt","ESD_TrueNegPion_Pt",1000,0,25);
 573                         fTrueList[iCut]->Add(fHistoTrueNegPionPt[iCut]);
 574
 575                         fHistoTrueNegPionFromNeutralMesonPt[iCut] = new TH1F("ESD_TrueNegPionFromNeutralMeson_Pt","ESD_TrueNegPionFromNeutralMeson_Pt",1000,0,25);
 576                         fTrueList[iCut]->Add(fHistoTrueNegPionFromNeutralMesonPt[iCut]);
 577                         fHistoTruePosPionFromNeutralMesonPt[iCut] = new TH1F("ESD_TruePosPionFromNeutralMeson_Pt","ESD_TruePosPionFromNeutralMeson_Pt",1000,0,25);
 578                         fTrueList[iCut]->Add(fHistoTruePosPionFromNeutralMesonPt[iCut]);
 579
 580                         fHistoTrueMotherPiPlPiMiPiZeroInvMassPt[iCut] = new TH2F("ESD_TrueMotherPiPlPiMiPiZero_InvMass_Pt","ESD_TrueMotherPiPlPiMiPiZero_InvMass_Pt",500,0.4,0.9,250,0,25);
 581                         fHistoTrueMotherPiPlPiMiPiZeroInvMassPt[iCut]->Sumw2();
 582                         fTrueList[iCut]->Add(fHistoTrueMotherPiPlPiMiPiZeroInvMassPt[iCut]);
 583
 584                         fHistoTrueMotherGammaGammaInvMassPt[iCut] = new TH2F("ESD_TrueMotherGG_InvMass_Pt","ESD_TrueMotherGG_InvMass_Pt",450,0.,0.45,250,0,25);
 585                         fHistoTrueMotherGammaGammaInvMassPt[iCut]->Sumw2();
 586                         fTrueList[iCut]->Add(fHistoTrueMotherGammaGammaInvMassPt[iCut]);
 587
 588 //                      if (fDoMesonQA){
 589                                 fHistoTruePionPionInvMassPt[iCut] = new TH2F("ESD_TruePiPlusPiNeg_InvMassPt","ESD_TruePiPlusPiNeg_InvMassPt",2000,0.,2.,200,0.,20.);
 590                                 fTrueList[iCut]->Add(fHistoTruePionPionInvMassPt[iCut]);
 591                                 fHistoTruePionPionFromNeutralMesonInvMassPt[iCut] = new TH2F("ESD_TruePiPlusPiNegFromNeutralMeson_InvMassPt","ESD_TruePiPlusPiNegFromNeutralMeson_InvMassPt",2000,0.,2.,200,0.,20.);
 592                                 fTrueList[iCut]->Add(fHistoTruePionPionFromNeutralMesonInvMassPt[iCut]);
 593 //                      }
 594                 }
 595         }
 596
 597
 598
 599         InitBack(); // Init Background Handler
 600
 601         fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
 602         if(!fV0Reader){printf("Error: No V0 Reader");return;} // GetV0Reader
 603
 604         if(fV0Reader)
 605                 if((AliConversionPhotonCuts*)fV0Reader->GetConversionCuts())
 606                         if(((AliConversionPhotonCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms())
 607                                 fOutputContainer->Add(((AliConversionPhotonCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms());
 608
 609
 610
 611         fPionSelector=(AliPrimaryPionSelector*)AliAnalysisManager::GetAnalysisManager()->GetTask("PionSelector");
 612         if(!fPionSelector){printf("Error: No PionSelector");return;} // GetV0Reader
 613
 614         if( fPionSelector ){
 615                 if ( ((AliPrimaryPionCuts*)fPionSelector->GetPrimaryPionCuts())->GetCutHistograms() ){
 616                         fOutputContainer->Add( ((AliPrimaryPionCuts*)fPionSelector->GetPrimaryPionCuts())->GetCutHistograms() );
 617                 }
 618         }
 619
 620         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 621                 if( fPionCutArray ){
 622                         if( ((AliPrimaryPionCuts*)fPionCutArray->At(iCut))->GetCutHistograms() ) {
 623                                 fCutFolder[iCut]->Add( ((AliPrimaryPionCuts*)fPionCutArray->At(iCut))->GetCutHistograms() );
 624                         }
 625                 }
 626                 if( fGammaCutArray ) {
 627                         if( ((AliConversionPhotonCuts*)fGammaCutArray->At(iCut))->GetCutHistograms() ) {
 628                                 fCutFolder[iCut]->Add( ((AliConversionPhotonCuts*)fGammaCutArray->At(iCut))->GetCutHistograms()  );
 629                         }
 630                 }
 631                 if( fNeutralPionMesonCutArray  ) {
 632                         if( ((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(iCut))->GetCutHistograms() ) {
 633                                 fCutFolder[iCut]->Add( ((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(iCut))->GetCutHistograms());
 634                         }
 635                 }
 636                 if( fMesonCutArray  ) {
 637                         if( ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutHistograms() ) {
 638                                 fCutFolder[iCut]->Add( ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutHistograms());
 639                         }
 640                 }
 641         }
 642
 643         PostData(1, fOutputContainer);
 644
 645 }
 646
 647 //______________________________________________________________________
 648 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::UserExec(Option_t *){
 649
 650         //
 651         // Execute analysis for current event
 652         //
 653
 654         fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
 655         if(!fV0Reader){printf("Error: No V0 Reader");return;} // GetV0Reader
 656
 657         Int_t eventQuality = ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEventQuality();
 658
 659         if(eventQuality == 2 || eventQuality == 3){// Event Not Accepted due to MC event missing or wrong trigger for V0ReaderV1
 660                 for(Int_t iCut = 0; iCut<fnCuts; iCut++){
 661                         fHistoNEvents[iCut]->Fill(eventQuality);
 662                 }
 663                 return;
 664         }
 665
 666         fPionSelector=(AliPrimaryPionSelector*)AliAnalysisManager::GetAnalysisManager()->GetTask("PionSelector");
 667         if(!fPionSelector){printf("Error: No PionSelector");return;} // GetV0Reader
 668
 669         if(fIsMC) fMCEvent     =  MCEvent();
 670         fESDEvent        = (AliESDEvent*)InputEvent();
 671         fReaderGammas    = fV0Reader->GetReconstructedGammas(); // Gammas from default Cut
 672         fSelectorNegPionIndex = fPionSelector->GetReconstructedNegPionIndex(); // Electrons from default Cut
 673         fSelectorPosPionIndex = fPionSelector->GetReconstructedPosPionIndex(); // Positrons from default Cut
 674
 675         fNumberOfESDTracks = fV0Reader->GetNumberOfPrimaryTracks();
 676         //AddTaskContainers(); //Add conatiner
 677
 678         for(Int_t iCut = 0; iCut<fnCuts; iCut++){
 679                 fiCut = iCut;
 680                 Int_t eventNotAccepted = ((AliConvEventCuts*)fEventCutArray->At(iCut))->IsEventAcceptedByCut(fV0Reader->GetEventCuts(),fInputEvent,fMCEvent,fIsHeavyIon);
 681
 682                 if(eventNotAccepted){
 683                         //                      cout << "event rejected due to wrong trigger: " <<eventNotAccepted << endl;
 684                         fHistoNEvents[iCut]->Fill(eventNotAccepted); // Check Centrality, PileUp, SDD and V0AND --> Not Accepted => eventQuality = 1
 685                         continue;
 686                 }
 687
 688                 if(eventQuality != 0){// Event Not Accepted
 689                         //                      cout << "event rejected due to: " <<eventQuality << endl;
 690                         fHistoNEvents[iCut]->Fill(eventQuality);
 691                         continue;
 692                 }
 693
 694                 fHistoNEvents[iCut]->Fill(eventQuality);
 695                 fHistoNGoodESDTracks[iCut]->Fill(fNumberOfESDTracks);
 696
 697                 if(fMCEvent){ // Process MC Particle
 698                         fMCStack = fMCEvent->Stack();
 699                         if(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetSignalRejection() != 0){
 700                                 ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetNotRejectedParticles(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetSignalRejection(),
 701                                                                                                                                                                                 ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetAcceptedHeader(),
 702                                                                                                                                                                                 fMCEvent);
 703                         }
 704                         ProcessMCParticles();
 705                 }
 706
 707                 ProcessPhotonCandidates(); // Process this cuts gammas
 708                 ProcessNeutralPionCandidatesPureConversions(); // Process neutral pion candidates
 709                 ProcessPionCandidates(); // Process this cuts gammas
 710
 711                 CalculateMesonCandidates();
 712                 CalculateBackground();
 713                 UpdateEventByEventData();
 714
 715                 fGoodGammas->Clear(); // delete this cuts good gammas
 716 //              if (fNeutralPionCandidates->GetEntries()>0)cout << "##################################"<<  fNeutralPionCandidates->GetEntries() << endl;
 717                 if(fNeutralPionCandidates->GetEntries()>0){
 718                         fNeutralPionCandidates->Clear();
 719                 }
 720                 fGoodVirtualParticles->Clear(); // delete this cuts good gammas
 721         }
 722
 723         fSelectorNegPionIndex.clear();
 724         fSelectorPosPionIndex.clear();
 725
 726         PostData( 1, fOutputContainer );
 727 }
 728
 729 Bool_t AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::Notify(){
 730         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 731                 if( !((AliConvEventCuts*)fEventCutArray->At(iCut))->GetDoEtaShift() ){
 732                         fProfileEtaShift[iCut]->Fill(0.,0.);
 733                         continue; // No Eta Shift requested, continue
 734                 }
 735                 if( ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift() == 0.0){ // Eta Shift requested but not set, get shift automatically
 736                         ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCorrectEtaShiftFromPeriod(fV0Reader->GetPeriodName());
 737                         ((AliConvEventCuts*)fEventCutArray->At(iCut))->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
 738                         ((AliPrimaryPionCuts*)fPionCutArray->At(iCut))->SetEtaShift( ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift() );
 739                         fProfileEtaShift[iCut]->Fill(0.,(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift()));
 740                         continue;
 741                 } else {
 742                         printf(" Eta t PiPlusPiMinus Gamma Task %s :: Eta Shift Manually Set to %f \n\n",
 743                         (((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber()).Data(),((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift());
 744                         ((AliConvEventCuts*)fEventCutArray->At(iCut))->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
 745                         ((AliPrimaryPionCuts*)fPionCutArray->At(iCut))->SetEtaShift( ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift() );
 746                         fProfileEtaShift[iCut]->Fill(0.,(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift()));
 747                 }
 748         }
 749         return kTRUE;
 750 }
 751
 752
 753 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::Terminate(const Option_t *){
 754 ///Grid
 755 }
 756
 757 //________________________________________________________________________
 758 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessPhotonCandidates(){
 759         Int_t nV0 = 0;
 760         TList *GoodGammasStepOne = new TList();
 761         TList *GoodGammasStepTwo = new TList();
 762         // Loop over Photon Candidates allocated by ReaderV1
 763
 764         for(Int_t i = 0; i < fReaderGammas->GetEntriesFast(); i++){
 765                 AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) fReaderGammas->At(i);
 766                 if(!PhotonCandidate) continue;
 767
 768                 fIsFromMBHeader = kTRUE;
 769
 770                 if( fMCEvent && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0 ){
 771                         Int_t isPosFromMBHeader
 772                                 = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack, fInputEvent);
 773                         if(isPosFromMBHeader == 0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
 774                         Int_t isNegFromMBHeader
 775                                 = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack,fInputEvent);
 776                         if(isNegFromMBHeader == 0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
 777                         if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE;
 778                 }
 779
 780                 if(!((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->PhotonIsSelected(PhotonCandidate,fESDEvent)) continue;
 781
 782                 if(!((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->UseElecSharingCut() &&
 783                         !((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->UseToCloseV0sCut()){ // if no post reader loop is required add to events good gammas
 784
 785                         fGoodGammas->Add(PhotonCandidate);
 786
 787                         if(fIsFromMBHeader){
 788                                 fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
 789                                 fHistoConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
 790                         }
 791
 792                         if(fMCEvent){
 793                                 ProcessTruePhotonCandidates(PhotonCandidate);
 794                         }
 795                 } else if(((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->UseElecSharingCut()){ // if Shared Electron cut is enabled, Fill array, add to step one
 796                         ((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->FillElectonLabelArray(PhotonCandidate,nV0);
 797                         nV0++;
 798                         GoodGammasStepOne->Add(PhotonCandidate);
 799                 } else if(!((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->UseElecSharingCut() &&
 800                                 ((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->UseToCloseV0sCut()){ // shared electron is disabled, step one not needed -> step two
 801                         GoodGammasStepTwo->Add(PhotonCandidate);
 802                 }
 803         }
 804
 805
 806         if(((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->UseElecSharingCut()){
 807                 for(Int_t i = 0;i<GoodGammasStepOne->GetEntries();i++){
 808                         AliAODConversionPhoton *PhotonCandidate= (AliAODConversionPhoton*) GoodGammasStepOne->At(i);
 809                         if(!PhotonCandidate) continue;
 810                         fIsFromMBHeader = kTRUE;
 811                         if(fMCEvent && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0){
 812                                 Int_t isPosFromMBHeader
 813                                 = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack,fInputEvent);
 814                                 Int_t isNegFromMBHeader
 815                                 = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack,fInputEvent);
 816                                 if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE;
 817                         }
 818                         if(!((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->RejectSharedElectronV0s(PhotonCandidate,i,GoodGammasStepOne->GetEntries())) continue;
 819                         if(!((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->UseToCloseV0sCut()){ // To Colse v0s cut diabled, step two not needed
 820                                 fGoodGammas->Add(PhotonCandidate);
 821                                 if(fIsFromMBHeader){
 822                                         fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
 823                                         fHistoConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
 824                                 }
 825                                 if(fMCEvent){
 826                                         ProcessTruePhotonCandidates(PhotonCandidate);
 827                                 }
 828                         }
 829                         else GoodGammasStepTwo->Add(PhotonCandidate); // Close v0s cut enabled -> add to list two
 830                 }
 831         }
 832         if(((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->UseToCloseV0sCut()){
 833                 for(Int_t i = 0;i<GoodGammasStepTwo->GetEntries();i++){
 834                         AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) GoodGammasStepTwo->At(i);
 835                         if(!PhotonCandidate) continue;
 836
 837                         if(fMCEvent && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0){
 838                                 Int_t isPosFromMBHeader
 839                                 = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack,fInputEvent);
 840                                 Int_t isNegFromMBHeader
 841                                 = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack,fInputEvent);
 842                                 if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE;
 843                         }
 844
 845                         if(!((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->RejectToCloseV0s(PhotonCandidate,GoodGammasStepTwo,i)) continue;
 846                         fGoodGammas->Add(PhotonCandidate); // Add gamma to current cut TList
 847
 848                         if(fIsFromMBHeader){
 849                                 fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt()); // Differences to old V0Reader in p_t due to conversion KF->TLorentzVector
 850                                 fHistoConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
 851                         }
 852
 853                         if(fMCEvent){
 854                                 ProcessTruePhotonCandidates(PhotonCandidate);
 855                         }
 856                 }
 857         }
 858
 859         delete GoodGammasStepOne;
 860         GoodGammasStepOne = 0x0;
 861         delete GoodGammasStepTwo;
 862         GoodGammasStepTwo = 0x0;
 863 }
 864
 865 //________________________________________________________________________
 866 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessTruePhotonCandidates(AliAODConversionPhoton *TruePhotonCandidate)
 867 {
 868         // Process True Photons
 869         AliStack *MCStack = fMCEvent->Stack();
 870         TParticle *posDaughter = TruePhotonCandidate->GetPositiveMCDaughter(MCStack);
 871         TParticle *negDaughter = TruePhotonCandidate->GetNegativeMCDaughter(MCStack);
 872
 873         if(posDaughter == NULL || negDaughter == NULL) return; // One particle does not exist
 874         if(posDaughter->GetMother(0) != negDaughter->GetMother(0)){  // Not Same Mother == Combinatorial Bck
 875                 return;
 876         }
 877
 878         else if (posDaughter->GetMother(0) == -1){
 879                 return;
 880         }
 881
 882         if(TMath::Abs(posDaughter->GetPdgCode())!=11 || TMath::Abs(negDaughter->GetPdgCode())!=11) return; //One Particle is not electron
 883         if(posDaughter->GetPdgCode()==negDaughter->GetPdgCode()) return; // Same Charge
 884         if(posDaughter->GetUniqueID() != 5 || negDaughter->GetUniqueID() !=5) return;// check if the daughters come from a conversion
 885
 886         TParticle *Photon = TruePhotonCandidate->GetMCParticle(MCStack);
 887         if(Photon->GetPdgCode() != 22) return; // Mother is no Photon
 888
 889         // True Photon
 890
 891         Int_t labelGamma = TruePhotonCandidate->GetMCParticleLabel(MCStack);
 892
 893         if( labelGamma < MCStack->GetNprimary() ){
 894                 if( fIsFromMBHeader ){
 895                         fHistoTrueConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
 896                         if (GammaIsNeutralMesonPiPlPiMiPiZeroDaughter(labelGamma)){
 897                                 fHistoTrueConvGammaFromNeutralMesonPt[fiCut]->Fill(TruePhotonCandidate->Pt());
 898                         }
 899                 }
 900         }
 901 }
 902
 903 //________________________________________________________________________
 904 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessNeutralPionCandidatesPureConversions(){
 905         // Conversion Gammas
 906         if(fGoodGammas->GetEntries()>1){
 907                 for(Int_t firstGammaIndex=0;firstGammaIndex<fGoodGammas->GetEntries()-1;firstGammaIndex++){
 908                         AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(firstGammaIndex));
 909                         if (gamma0==NULL) continue;
 910                         for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fGoodGammas->GetEntries();secondGammaIndex++){
 911                                 AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(secondGammaIndex));
 912                                 //Check for same Electron ID
 913                                 if (gamma1==NULL) continue;
 914                                 if(gamma0->GetTrackLabelPositive() == gamma1->GetTrackLabelPositive() ||
 915                                 gamma0->GetTrackLabelNegative() == gamma1->GetTrackLabelNegative() ||
 916                                 gamma0->GetTrackLabelNegative() == gamma1->GetTrackLabelPositive() ||
 917                                 gamma0->GetTrackLabelPositive() == gamma1->GetTrackLabelNegative() ) continue;
 918
 919                                 AliAODConversionMother *pi0cand = new AliAODConversionMother(gamma0,gamma1);
 920                                 pi0cand->SetLabels(firstGammaIndex,secondGammaIndex);
 921
 922                                 pi0cand->CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
 923                                 if((((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(fiCut))->MesonIsSelected(pi0cand,kTRUE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift()))){
 924                                         fHistoGammaGammaInvMassPt[fiCut]->Fill(pi0cand->M(),pi0cand->Pt());
 925                                         if (pi0cand->M() > ((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(fiCut))->GetSelectionLow() && pi0cand->M() < ((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(fiCut))->GetSelectionHigh()){
 926                                                 fNeutralPionCandidates->Add(pi0cand);
 927 //                                              cout << "Pi0 candidate " << pi0cand->M() << "\t" << pi0cand->Pt() << endl;
 928                                         }
 929
 930                                         if(fIsMC){
 931                                                 if(fInputEvent->IsA()==AliESDEvent::Class())
 932                                                         ProcessTrueNeutralPionCandidatesPureConversions(pi0cand,gamma0,gamma1);
 933                                                 if(fInputEvent->IsA()==AliAODEvent::Class())
 934                                                         ProcessTrueNeutralPionCandidatesPureConversionsAOD(pi0cand,gamma0,gamma1);
 935                                         }
 936                                 }
 937 //                              delete pi0cand;
 938 //                              pi0cand=0x0;
 939                         }
 940                 }
 941         }
 942 }
 943
 944
 945 //______________________________________________________________________
 946 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessTrueNeutralPionCandidatesPureConversions(AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1)
 947 {
 948         // Process True Mesons
 949         AliStack *MCStack = fMCEvent->Stack();
 950         if(TrueGammaCandidate0->GetV0Index()<fInputEvent->GetNumberOfV0s()){
 951                 Bool_t isTruePi0 = kFALSE;
 952                 Bool_t isTruePi0Dalitz = kFALSE;
 953                 Bool_t gamma0DalitzCand = kFALSE;
 954                 Bool_t gamma1DalitzCand = kFALSE;
 955                 Int_t gamma0MCLabel = TrueGammaCandidate0->GetMCParticleLabel(MCStack);
 956                 Int_t gamma0MotherLabel = -1;
 957                 Int_t motherRealLabel = -1;
 958                 if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
 959                         // Daughters Gamma 0
 960                         TParticle * negativeMC = (TParticle*)TrueGammaCandidate0->GetNegativeMCDaughter(MCStack);
 961                         TParticle * positiveMC = (TParticle*)TrueGammaCandidate0->GetPositiveMCDaughter(MCStack);
 962                         TParticle * gammaMC0 = (TParticle*)MCStack->Particle(gamma0MCLabel);
 963                         if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
 964                                 if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
 965                                         if(gammaMC0->GetPdgCode() == 22){ // ... with Gamma Mother
 966                                                 gamma0MotherLabel=gammaMC0->GetFirstMother();
 967                                                 motherRealLabel=gammaMC0->GetFirstMother();
 968                                         }
 969                                 }
 970                                 if(gammaMC0->GetPdgCode() ==111){ // Dalitz candidate
 971                                         gamma0DalitzCand = kTRUE;
 972                                         gamma0MotherLabel=-111;
 973                                 }
 974                         }
 975                 }
 976                 if(TrueGammaCandidate1->GetV0Index()<fInputEvent->GetNumberOfV0s()){
 977                         Int_t gamma1MCLabel = TrueGammaCandidate1->GetMCParticleLabel(MCStack);
 978                         Int_t gamma1MotherLabel = -1;
 979                         if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
 980                                 // Daughters Gamma 1
 981                                 TParticle * negativeMC = (TParticle*)TrueGammaCandidate1->GetNegativeMCDaughter(MCStack);
 982                                 TParticle * positiveMC = (TParticle*)TrueGammaCandidate1->GetPositiveMCDaughter(MCStack);
 983                                 TParticle * gammaMC1 = (TParticle*)MCStack->Particle(gamma1MCLabel);
 984                                 if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
 985                                         if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
 986                                                 if(gammaMC1->GetPdgCode() == 22){ // ... with Gamma Mother
 987                                                         gamma1MotherLabel=gammaMC1->GetFirstMother();
 988                                                 }
 989                                         }
 990                                         if(gammaMC1->GetPdgCode() ==111 ){ // Dalitz candidate
 991                                                 gamma1DalitzCand = kTRUE;
 992                                                 gamma1MotherLabel=-111;
 993                                         }
 994                                 }
 995                         }
 996                         if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
 997                                 if(((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetPdgCode() == 111){
 998                                         isTruePi0=kTRUE;
 999                                 }
1000                         }
1001
1002                         //Identify Dalitz candidate
1003                         if (gamma1DalitzCand || gamma0DalitzCand){
1004                                 if (gamma0DalitzCand && gamma0MCLabel >=0 && gamma0MCLabel==gamma1MotherLabel){
1005                                         if (gamma0MotherLabel == -111) isTruePi0Dalitz = kTRUE;
1006                                 }
1007                                 if (gamma1DalitzCand && gamma1MCLabel >=0 && gamma1MCLabel==gamma0MotherLabel){
1008                                         if (gamma1MotherLabel == -111) isTruePi0Dalitz = kTRUE;
1009                                 }
1010                         }
1011
1012
1013                         if(isTruePi0 || isTruePi0Dalitz){// True Pion
1014                                 Pi0Candidate->SetTrueMesonValue(1);
1015                                 Pi0Candidate->SetMCLabel(motherRealLabel);
1016                                 fHistoTrueMotherGammaGammaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1017                         }
1018                 }
1019         }
1020 }
1021 //______________________________________________________________________
1022 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessTrueNeutralPionCandidatesPureConversionsAOD(AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1)
1023 {
1024
1025         // Process True Mesons
1026         TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1027         Bool_t isTruePi0 = kFALSE;
1028         Bool_t isTruePi0Dalitz = kFALSE;
1029         Bool_t gamma0DalitzCand = kFALSE;
1030         Bool_t gamma1DalitzCand = kFALSE;
1031         Int_t motherRealLabel = -1;
1032
1033         if (AODMCTrackArray!=NULL && TrueGammaCandidate0 != NULL){
1034                 AliAODMCParticle *positiveMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate0->GetMCLabelPositive()));
1035                 AliAODMCParticle *negativeMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate0->GetMCLabelNegative()));
1036
1037                 Int_t gamma0MCLabel = -1;
1038                 Int_t gamma0MotherLabel = -1;
1039                 if(!positiveMC||!negativeMC)
1040                         return;
1041
1042                 if(positiveMC->GetMother()>-1&&(negativeMC->GetMother() == positiveMC->GetMother())){
1043                         gamma0MCLabel = positiveMC->GetMother();
1044                 }
1045
1046                 if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1047                         // Daughters Gamma 0
1048                         AliAODMCParticle * gammaMC0 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma0MCLabel));
1049                         if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1050                                 if(((positiveMC->GetMCProcessCode())) == 5 && ((negativeMC->GetMCProcessCode())) == 5){ // ... From Conversion ...
1051                                         if(gammaMC0->GetPdgCode() == 22){ // ... with Gamma Mother
1052                                                 gamma0MotherLabel=gammaMC0->GetMother();
1053                                                 motherRealLabel=gammaMC0->GetMother();
1054                                         }
1055                                 }
1056                                 if(gammaMC0->GetPdgCode() ==111){ // Dalitz candidate
1057                                         gamma0DalitzCand = kTRUE;
1058                                         gamma0MotherLabel=-111;
1059                                 }
1060                         }
1061                 }
1062                 positiveMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate1->GetMCLabelPositive()));
1063                 negativeMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate1->GetMCLabelNegative()));
1064
1065                 Int_t gamma1MCLabel = -1;
1066                 Int_t gamma1MotherLabel = -1;
1067                 if(!positiveMC||!negativeMC)
1068                         return;
1069
1070                 if(positiveMC->GetMother()>-1&&(negativeMC->GetMother() == positiveMC->GetMother())){
1071                         gamma1MCLabel = positiveMC->GetMother();
1072                 }
1073                 if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1074                         // Daughters Gamma 1
1075                         AliAODMCParticle * gammaMC1 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MCLabel));
1076                         if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1077                                 if(((positiveMC->GetMCProcessCode())) == 5 && ((negativeMC->GetMCProcessCode())) == 5){ // ... From Conversion ...
1078                                         if(gammaMC1->GetPdgCode() == 22){ // ... with Gamma Mother
1079                                         gamma1MotherLabel=gammaMC1->GetMother();
1080                                         }
1081                                 }
1082                                 if(gammaMC1->GetPdgCode() ==111 ){ // Dalitz candidate
1083                                                 gamma1DalitzCand = kTRUE;
1084                                                 gamma1MotherLabel=-111;
1085                                 }
1086                         }
1087                 }
1088                 if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
1089                         if(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->GetPdgCode() == 111){
1090                                 isTruePi0=kTRUE;
1091                         }
1092                 }
1093
1094                 //Identify Dalitz candidate
1095                 if (gamma1DalitzCand || gamma0DalitzCand){
1096                         if (gamma0DalitzCand && gamma0MCLabel >=0 && gamma0MCLabel==gamma1MotherLabel){
1097                                 if (gamma0MotherLabel == -111) isTruePi0Dalitz = kTRUE;
1098                         }
1099                         if (gamma1DalitzCand && gamma1MCLabel >=0 && gamma1MCLabel==gamma0MotherLabel){
1100                                 if (gamma1MotherLabel == -111) isTruePi0Dalitz = kTRUE;
1101                         }
1102                 }
1103
1104                 if(isTruePi0 || isTruePi0Dalitz){// True Pion
1105                         Pi0Candidate->SetTrueMesonValue(1);
1106                         Pi0Candidate->SetMCLabel(motherRealLabel);
1107                         fHistoTrueMotherGammaGammaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1108                 }
1109         }
1110         return;
1111 }
1112
1113
1114
1115 //________________________________________________________________________
1116 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessPionCandidates(){
1117
1118         Double_t magField = fInputEvent->GetMagneticField();
1119         if( magField  < 0.0 ){
1120                 magField =  1.0;
1121         } else {
1122                 magField =  -1.0;
1123         }
1124
1125         vector<Int_t> lGoodNegPionIndexPrev(0);
1126         vector<Int_t> lGoodPosPionIndexPrev(0);
1127
1128         for(UInt_t i = 0; i < fSelectorNegPionIndex.size(); i++){
1129                 AliESDtrack* negPionCandidate = fESDEvent->GetTrack(fSelectorNegPionIndex[i]);
1130                 if(! ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->PionIsSelected(negPionCandidate) ) continue;
1131                 lGoodNegPionIndexPrev.push_back(   fSelectorNegPionIndex[i] );
1132                 fHistoNegPionPt[fiCut]->Fill(negPionCandidate->Pt());
1133                 fHistoNegPionPhi[fiCut]->Fill(negPionCandidate->Phi());
1134                 if( fMCEvent ) {
1135                         Int_t labelNegPion = TMath::Abs( negPionCandidate->GetLabel() );
1136                         if( labelNegPion < fMCStack->GetNtrack() ){
1137                                 TParticle* negPion = fMCStack->Particle(labelNegPion);
1138                                 if( negPion->GetPdgCode() ==  -211 ){
1139                                         if( labelNegPion < fMCStack->GetNprimary() ){
1140                                                 fHistoTrueNegPionPt[fiCut]->Fill(negPionCandidate->Pt());    //primary negPion
1141                                         }
1142                                         if( IsEtaPiPlPiMiPiZeroDaughter(labelNegPion) || IsOmegaPiPlPiMiPiZeroDaughter(labelNegPion) ) {
1143                                                 if( labelNegPion < fMCStack->GetNprimary() ) {
1144                                                         fHistoTrueNegPionFromNeutralMesonPt[fiCut]->Fill(negPionCandidate->Pt());
1145                                                 }
1146                                         }
1147                                 }
1148                         }
1149                 }
1150         }
1151
1152         for(UInt_t i = 0; i < fSelectorPosPionIndex.size(); i++){
1153                 AliESDtrack* posPionCandidate = fESDEvent->GetTrack( fSelectorPosPionIndex[i] );
1154                 if(! ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->PionIsSelected(posPionCandidate) ) continue;
1155                 lGoodPosPionIndexPrev.push_back(   fSelectorPosPionIndex[i]  );
1156                 fHistoPosPionPt[fiCut]->Fill( posPionCandidate->Pt() );
1157                 fHistoPosPionPhi[fiCut]->Fill( posPionCandidate->Phi() );
1158
1159                 if( fMCEvent ) {
1160                         Int_t labelPosPion = TMath::Abs( posPionCandidate->GetLabel() );
1161                         if( labelPosPion < fMCStack->GetNtrack() ) {
1162                                 TParticle* posPion = fMCStack->Particle(labelPosPion);
1163                                 if( posPion->GetPdgCode() ==  211 ){
1164                                         if( labelPosPion < fMCStack->GetNprimary() ){
1165                                                 fHistoTruePosPionPt[fiCut]->Fill(posPionCandidate->Pt());
1166                                         }
1167                                         if( IsEtaPiPlPiMiPiZeroDaughter(labelPosPion) || IsOmegaPiPlPiMiPiZeroDaughter(labelPosPion) ) {
1168                                                 if( labelPosPion < fMCStack->GetNprimary() ){
1169                                                         fHistoTruePosPionFromNeutralMesonPt[fiCut]->Fill(posPionCandidate->Pt());
1170                                                 }
1171                                         }
1172                                 }
1173                         }
1174                 }
1175         }
1176
1177
1178         for(UInt_t i = 0; i < lGoodNegPionIndexPrev.size(); i++){
1179
1180                 AliESDtrack *negPionCandidate = fESDEvent->GetTrack(lGoodNegPionIndexPrev[i]);
1181                 AliKFParticle negPionCandidateKF( *negPionCandidate->GetConstrainedParam(), 211 );
1182
1183                 for(UInt_t j = 0; j < lGoodPosPionIndexPrev.size(); j++){
1184
1185                         AliESDtrack *posPionCandidate = fESDEvent->GetTrack(lGoodPosPionIndexPrev[j]);
1186                         AliKFParticle posPionCandidateKF( *posPionCandidate->GetConstrainedParam(), 211 );
1187
1188                         AliKFConversionPhoton* virtualPhoton = NULL;
1189                         virtualPhoton = new AliKFConversionPhoton(negPionCandidateKF,posPionCandidateKF);
1190                         AliKFVertex primaryVertexImproved(*fInputEvent->GetPrimaryVertex());
1191 //                      primaryVertexImproved+=*virtualPhoton;
1192                         virtualPhoton->SetProductionVertex(primaryVertexImproved);
1193                         virtualPhoton->SetTrackLabels( lGoodPosPionIndexPrev[j], lGoodNegPionIndexPrev[i]);
1194
1195
1196                         if( fMCEvent ) {
1197                                 Int_t labeln=TMath::Abs(negPionCandidate->GetLabel());
1198                                 Int_t labelp=TMath::Abs(posPionCandidate->GetLabel());
1199                                 TParticle *fNegativeMCParticle = fMCStack->Particle(labeln);
1200                                 TParticle *fPositiveMCParticle = fMCStack->Particle(labelp);
1201
1202                                 if( fPositiveMCParticle && fNegativeMCParticle) {
1203                                         virtualPhoton->SetMCLabelPositive(labelp);
1204                                         virtualPhoton->SetMCLabelNegative(labeln);
1205                                 }
1206                         }
1207
1208                         AliAODConversionPhoton *vParticle = new AliAODConversionPhoton(virtualPhoton); //To Apply PsiPairCut
1209                         if (((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->DoMassCut()){
1210                                 if (vParticle->GetMass() < ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->GetMassCut()){
1211                                         fGoodVirtualParticles->Add(  vParticle );
1212                                 }
1213                         } else {
1214                                 fGoodVirtualParticles->Add(  vParticle );
1215                         }
1216                         delete virtualPhoton;
1217                         virtualPhoton=NULL;
1218
1219                 }
1220         }
1221 }
1222
1223 //_____________________________________________________________________________
1224 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessMCParticles(){
1225
1226         // Loop over all primary MC particle
1227
1228         for(Int_t i = 0; i < fMCStack->GetNprimary(); i++) {
1229
1230                 TParticle* particle = (TParticle *)fMCStack->Particle(i);
1231                 if (!particle) continue;
1232
1233                 Int_t isMCFromMBHeader = -1;
1234                 if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0){
1235                         isMCFromMBHeader
1236                                 = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent);
1237                         if(isMCFromMBHeader == 0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1238                 }
1239
1240                 if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack,fInputEvent)){
1241
1242                         // find MC photons
1243                         if(((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->PhotonIsSelectedMC(particle,fMCStack,kFALSE)){
1244                                 fHistoMCAllGammaPt[fiCut]->Fill(particle->Pt()); // All MC Gamma
1245                                 if(particle->GetMother(0) >-1){
1246                                         if (fMCStack->Particle(particle->GetMother(0))->GetPdgCode() ==111){
1247                                                 if (fMCStack->Particle(particle->GetMother(0))->GetMother(0) > -1){
1248                                                         if ( fMCStack->Particle((fMCStack->Particle(particle->GetMother(0)))->GetMother(0))->GetPdgCode() == 221 ||
1249                                                                  fMCStack->Particle((fMCStack->Particle(particle->GetMother(0)))->GetMother(0))->GetPdgCode() == 223 ){
1250                                                                 if ( fMCStack->Particle(particle->GetMother(0))->GetNDaughters()==3 )
1251                                                                         fHistoMCGammaFromNeutralMesonPt[fiCut]->Fill(particle->Pt()); // All photons from eta or omega via pi0
1252                                                         }
1253                                                 }
1254                                         }
1255                                 }
1256                         }
1257
1258                         if(((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->PhotonIsSelectedMC(particle,fMCStack,kTRUE)){
1259                                 fHistoMCConvGammaPt[fiCut]->Fill(particle->Pt());
1260                         } // Converted MC Gamma
1261
1262                         if(((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->PionIsSelectedMC(i,fMCStack)){
1263                                 if( particle->GetPdgCode() == 211){
1264                                         fHistoMCAllPosPionsPt[fiCut]->Fill(particle->Pt()); // All pos pions
1265                                         if(particle->GetMother(0) >-1){
1266                                                 if (fMCStack->Particle(particle->GetMother(0))->GetPdgCode() ==221 || fMCStack->Particle(particle->GetMother(0))->GetPdgCode() ==223)
1267                                                         fHistoMCPosPionsFromNeutralMesonPt[fiCut]->Fill(particle->Pt()); // All pos from eta or omega
1268                                         }
1269                                 }
1270                                 if( particle->GetPdgCode() == -211){
1271                                         fHistoMCAllNegPionsPt[fiCut]->Fill(particle->Pt()); // All neg pions
1272                                         if(particle->GetMother(0) >-1){
1273                                                 if (fMCStack->Particle(particle->GetMother(0))->GetPdgCode() ==221 || fMCStack->Particle(particle->GetMother(0))->GetPdgCode() ==223 )
1274                                                         fHistoMCNegPionsFromNeutralMesonPt[fiCut]->Fill(particle->Pt()); // All pos from eta or omega
1275                                         }
1276                                 }
1277                         }
1278
1279
1280                         // \eta -> pi+ pi- \gamma
1281                         Int_t labelNeutPion = -1;
1282                         Int_t labelNegPion = -1;
1283                         Int_t labelPosPion = -1;
1284
1285                         if( ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelectedMCPiPlPiMiPiZero(particle,fMCStack,labelNegPion,labelPosPion,labelNeutPion,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift())){
1286                                 Float_t weighted= 1;
1287                                 if( ((AliPrimaryPionCuts*) fPionCutArray->At(fiCut))->DoWeights() ) {
1288                                         if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack,fInputEvent)){
1289                                                 if (particle->Pt()>0.005){
1290                                                         weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, fMCStack,fInputEvent);
1291                                                 }
1292                                         }
1293                                 }
1294                                 if(particle->GetPdgCode() == 221)fHistoMCEtaPiPlPiMiPiZeroPt[fiCut]->Fill(particle->Pt(), weighted);                                            // All MC Eta in respective decay channel
1295                                 if(particle->GetPdgCode() == 223)fHistoMCOmegaPiPlPiMiPiZeroPt[fiCut]->Fill(particle->Pt(), weighted);                                          // All MC Omega in respective decay channel
1296
1297                                 TParticle *neutPion    = fMCStack->Particle(labelNeutPion);
1298                                 TParticle *gamma1 = fMCStack->Particle(neutPion->GetDaughter(0));
1299                                 TParticle *gamma2 = fMCStack->Particle(neutPion->GetDaughter(1));
1300                                 if(
1301                                         ((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->PhotonIsSelectedMC(gamma1,fMCStack,kFALSE) &&                                    // test first daugther of pi0
1302                                         ((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->PhotonIsSelectedMC(gamma2,fMCStack,kFALSE) &&                                    // test second daughter of pi0
1303                                         ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->PionIsSelectedMC(labelNegPion,fMCStack) &&                                                             // test negative pion
1304                                         ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->PionIsSelectedMC(labelPosPion,fMCStack)                                                                // test positive pion
1305                                    ) {
1306                                                 if(particle->GetPdgCode() == 221) fHistoMCEtaPiPlPiMiPiZeroInAccPt[fiCut]->Fill(particle->Pt(), weighted );             // MC Eta pi+ pi- pi0 with gamma's and e+e- in acc
1307                                                 if(particle->GetPdgCode() == 223) fHistoMCOmegaPiPlPiMiPiZeroInAccPt[fiCut]->Fill(particle->Pt(), weighted );           // MC Omega pi+ pi- pi0 with gamma's and e+e- in acc
1308                                 }
1309                         }
1310                 }
1311         }
1312 }
1313
1314
1315 //________________________________________________________________________
1316 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::CalculateMesonCandidates(){
1317
1318         // Conversion Gammas
1319         if( fNeutralPionCandidates->GetEntries() > 0 && fGoodVirtualParticles->GetEntries() > 0 ){
1320
1321                 vector<Bool_t> lGoodVirtualParticle(fGoodVirtualParticles->GetEntries(), kFALSE);
1322
1323                 for(Int_t mesonIndex=0; mesonIndex<fNeutralPionCandidates->GetEntries(); mesonIndex++){
1324                         AliAODConversionMother *neutralPion=dynamic_cast<AliAODConversionMother*>(fNeutralPionCandidates->At(mesonIndex));
1325                         if (neutralPion==NULL) continue;
1326                         for(Int_t virtualParticleIndex=0;virtualParticleIndex<fGoodVirtualParticles->GetEntries();virtualParticleIndex++){
1327
1328                                 AliAODConversionPhoton *vParticle=dynamic_cast<AliAODConversionPhoton*>(fGoodVirtualParticles->At(virtualParticleIndex));
1329                                 if (vParticle==NULL) continue;
1330                                 //Check for same Electron ID
1331
1332                                 AliAODConversionMother *mesoncand = new AliAODConversionMother(neutralPion,vParticle);
1333                                 mesoncand->SetLabels(mesonIndex,virtualParticleIndex);
1334
1335                                 if( ( ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(mesoncand,kTRUE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift())) ){
1336
1337 //                                      cout<< "Meson Accepted "<<endl;
1338
1339                                         Int_t zbin= fBGHandler[fiCut]->GetZBinIndex(fESDEvent->GetPrimaryVertex()->GetZ());
1340                                         Int_t mbin = 0;
1341                                         if( ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
1342                                                 mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fNumberOfESDTracks);
1343                                         } else {
1344                                                 mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGoodGammas->GetEntries());
1345                                         }
1346
1347                                         AliESDtrack *posPionVParticle = 0;
1348                                         AliESDtrack *negPionVParticle = 0;
1349
1350                                         Double_t clsToFPos = -1.0;
1351                                         Double_t clsToFNeg = -1.0;
1352
1353                                         Float_t dcaToVertexXYPos = -1.0;
1354                                         Float_t dcaToVertexZPos  = -1.0;
1355                                         Float_t dcaToVertexXYNeg = -1.0;
1356                                         Float_t dcaToVertexZNeg  = -1.0;
1357
1358
1359                                         if ( fDoMesonQA ) {
1360
1361                                                 fHistoPionPionInvMassPt[fiCut]->Fill( vParticle->GetMass(),vParticle->Pt());
1362
1363                                                 posPionVParticle = fESDEvent->GetTrack( vParticle->GetTrackLabelPositive() );
1364                                                 negPionVParticle = fESDEvent->GetTrack( vParticle->GetTrackLabelNegative() );
1365                                                 clsToFPos = ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->GetNFindableClustersTPC(posPionVParticle);
1366                                                 clsToFNeg = ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->GetNFindableClustersTPC(negPionVParticle);
1367
1368                                                 Float_t bPos[2];
1369                                                 Float_t bCovPos[3];
1370                                                 posPionVParticle->GetImpactParameters(bPos,bCovPos);
1371                                                 if (bCovPos[0]<=0 || bCovPos[2]<=0) {
1372                                                         AliDebug(1, "Estimated b resolution lower or equal zero!");
1373                                                         bCovPos[0]=0; bCovPos[2]=0;
1374                                                 }
1375
1376                                                 Float_t bNeg[2];
1377                                                 Float_t bCovNeg[3];
1378                                                 posPionVParticle->GetImpactParameters(bNeg,bCovNeg);
1379                                                 if (bCovNeg[0]<=0 || bCovNeg[2]<=0) {
1380                                                         AliDebug(1, "Estimated b resolution lower or equal zero!");
1381                                                         bCovNeg[0]=0; bCovNeg[2]=0;
1382                                                 }
1383
1384                                                 dcaToVertexXYPos = bPos[0];
1385                                                 dcaToVertexZPos  = bPos[1];
1386                                                 dcaToVertexXYNeg = bNeg[0];
1387                                                 dcaToVertexZNeg  = bNeg[1];
1388                                         }
1389
1390                                         fHistoMotherInvMassPt[fiCut]->Fill(mesoncand->M(),mesoncand->Pt());
1391                                         Double_t sparesFill[4] = {mesoncand->M(),mesoncand->Pt(),(Double_t)zbin,(Double_t)mbin};
1392                                         fTHnSparseMotherInvMassPtZM[fiCut]->Fill(sparesFill,1);
1393
1394                                         if ( fDoMesonQA ) {
1395                                                 if( lGoodVirtualParticle[virtualParticleIndex] == kFALSE ) {
1396
1397                                                         fHistoNegPionEta[fiCut]->Fill( negPionVParticle->Eta() );
1398                                                         fHistoPosPionEta[fiCut]->Fill( posPionVParticle->Eta() );
1399
1400                                                         fHistoNegPionClsTPC[fiCut]->Fill(clsToFNeg,negPionVParticle->Pt());
1401                                                         fHistoPosPionClsTPC[fiCut]->Fill(clsToFPos,posPionVParticle->Pt());
1402
1403                                                         fHistoPionDCAxy[fiCut]->Fill(  dcaToVertexXYNeg, negPionVParticle->Pt() );
1404                                                         fHistoPionDCAz[fiCut]->Fill(   dcaToVertexZNeg,  negPionVParticle->Pt() );
1405                                                         fHistoPionDCAxy[fiCut]->Fill(  dcaToVertexXYPos, posPionVParticle->Pt() );
1406                                                         fHistoPionDCAz[fiCut]->Fill(   dcaToVertexZPos,  posPionVParticle->Pt() );
1407
1408                                                         fHistoPionTPCdEdxNSigma[fiCut]->Fill( posPionVParticle->P(),((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->GetPIDResponse()->NumberOfSigmasTPC(posPionVParticle, AliPID::kPion) );
1409                                                         fHistoPionTPCdEdxNSigma[fiCut]->Fill( negPionVParticle->P(),((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->GetPIDResponse()->NumberOfSigmasTPC(negPionVParticle, AliPID::kPion) );
1410
1411                                                         fHistoPionTPCdEdx[fiCut]->Fill( posPionVParticle->P(), TMath::Abs(posPionVParticle->GetTPCsignal()));
1412                                                         fHistoPionTPCdEdx[fiCut]->Fill( negPionVParticle->P(), TMath::Abs(negPionVParticle->GetTPCsignal()));
1413
1414                                                         lGoodVirtualParticle[virtualParticleIndex] = kTRUE;
1415
1416                                                 }
1417                                         }
1418
1419
1420                                         if(fMCEvent){
1421                                                 ProcessTrueMesonCandidates(mesoncand,neutralPion,vParticle);
1422                                         }
1423                                 }
1424                                 delete mesoncand;
1425                                 mesoncand=0x0;
1426                         }
1427                 }
1428         }
1429 }
1430
1431 //________________________________________________________________________
1432 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::CalculateBackground(){
1433
1434         Int_t zbin= fBGHandler[fiCut]->GetZBinIndex(fESDEvent->GetPrimaryVertex()->GetZ());
1435         Int_t mbin = 0;
1436
1437
1438         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
1439                 mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fNumberOfESDTracks);
1440         } else {
1441                 mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGoodGammas->GetEntries());
1442         }
1443
1444         Int_t method = 1;
1445         AliGammaConversionAODBGHandler::GammaConversionVertex *bgEventVertex = NULL;
1446         if( ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity() ) {
1447                 for(Int_t nEventsInBG=0;nEventsInBG<fBGHandler[fiCut]->GetNBGEvents();nEventsInBG++){
1448                         AliGammaConversionMotherAODVector *previousEventMesons = fBGHandler[fiCut]->GetBGGoodMesons(zbin,mbin,nEventsInBG);
1449                         if(fMoveParticleAccordingToVertex == kTRUE && method == 1){
1450                                 bgEventVertex = fBGHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG);
1451                         }
1452
1453                         for(Int_t iCurrent=0;iCurrent<fGoodVirtualParticles->GetEntries();iCurrent++){
1454                                 AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGoodVirtualParticles->At(iCurrent));
1455
1456                                 for(UInt_t iPrevious=0;iPrevious<previousEventMesons->size();iPrevious++){
1457                                         AliAODConversionMother previousGoodMeson = (AliAODConversionMother)(*(previousEventMesons->at(iPrevious)));
1458
1459                                         if(fMoveParticleAccordingToVertex == kTRUE && method == 1 ){
1460                                                 MoveParticleAccordingToVertex(&previousGoodMeson,bgEventVertex);
1461                                         }
1462
1463                                         AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&previousGoodMeson,&currentEventGoodV0);
1464
1465
1466                                         if( ( ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE, ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift()))){
1467                                                 fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
1468                                                 Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin};
1469                                                 fTHnSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
1470                                         }
1471                                         delete backgroundCandidate;
1472                                         backgroundCandidate = 0x0;
1473                                 }
1474                         }
1475                 }
1476         } else {
1477                 for(Int_t nEventsInBG=0;nEventsInBG <fBGHandler[fiCut]->GetNBGEvents();nEventsInBG++){
1478                         AliGammaConversionMotherAODVector *previousEventMesons = fBGHandler[fiCut]->GetBGGoodMesons(zbin,mbin,nEventsInBG);
1479                         if(previousEventMesons){
1480                                 if(fMoveParticleAccordingToVertex == kTRUE && method == 1){
1481                                         bgEventVertex = fBGHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG);
1482                                 }
1483                                 for(Int_t iCurrent=0;iCurrent<fGoodVirtualParticles->GetEntries();iCurrent++){
1484
1485                                         AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGoodVirtualParticles->At(iCurrent));
1486
1487                                         for(UInt_t iPrevious=0;iPrevious<previousEventMesons->size();iPrevious++){
1488
1489                                                 AliAODConversionMother previousGoodMeson = (AliAODConversionMother)(*(previousEventMesons->at(iPrevious)));
1490
1491                                                 if(fMoveParticleAccordingToVertex == kTRUE && method ==1){
1492                                                         MoveParticleAccordingToVertex(&previousGoodMeson,bgEventVertex);
1493                                                 }
1494
1495                                                 AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&previousGoodMeson,&currentEventGoodV0);
1496
1497                                                 if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift()))){
1498                                                         fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
1499                                                         Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin};
1500                                                         fTHnSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
1501                                                 }
1502                                                 delete backgroundCandidate;
1503                                                 backgroundCandidate = 0x0;
1504                                         }
1505                                 }
1506                         }
1507                 }
1508         }
1509 }
1510
1511 //______________________________________________________________________
1512 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessTrueMesonCandidates(AliAODConversionMother *mesoncand, AliAODConversionMother *TrueNeutralPionCandidate, AliAODConversionPhoton *TrueVirtualParticleCandidate){
1513
1514         // Process True Mesons
1515         AliStack *MCStack = fMCEvent->Stack();
1516
1517         Bool_t isTrueEta = kFALSE;
1518         Bool_t isTrueOmega = kFALSE;
1519         Int_t trueMesonFlag = TrueNeutralPionCandidate->GetTrueMesonValue();
1520         Int_t pi0MCLabel= TrueNeutralPionCandidate->GetMCLabel();
1521
1522
1523         if ( !(trueMesonFlag == 1 && pi0MCLabel != -1)) return;
1524 //      cout << trueMesonFlag << "\t" << pi0MCLabel << endl;
1525
1526
1527         Int_t virtualParticleMCLabel = TrueVirtualParticleCandidate->GetMCParticleLabel(MCStack);
1528         Int_t virtualParticleMotherLabel = -1;
1529
1530         Bool_t isPiPiDecay = kFALSE;
1531
1532 //      if (fDoMesonQA){
1533                 TParticle * negativeMC = (TParticle*)TrueVirtualParticleCandidate->GetNegativeMCDaughter(MCStack);
1534                 TParticle * positiveMC = (TParticle*)TrueVirtualParticleCandidate->GetPositiveMCDaughter(MCStack);
1535                 if(TMath::Abs(negativeMC->GetPdgCode())==211 && TMath::Abs(positiveMC->GetPdgCode())==211){  // Pions ...
1536                         fHistoTruePionPionInvMassPt[fiCut]->Fill(TrueVirtualParticleCandidate->GetMass(),TrueVirtualParticleCandidate->Pt());
1537                 }
1538 //      }
1539
1540         if(virtualParticleMCLabel != -1){ // if virtualParticleMCLabel==-1 particles don't have same mother
1541 //              TParticle * negativeMC = (TParticle*)TrueVirtualParticleCandidate->GetNegativeMCDaughter(MCStack);
1542 //              TParticle * positiveMC = (TParticle*)TrueVirtualParticleCandidate->GetPositiveMCDaughter(MCStack);
1543 //                      TParticle * virtualParticleMotherMC = (TParticle*)MCStack->Particle(virtualParticleMCLabel);
1544 //                      cout << "pdg code same mother - " << virtualParticleMotherMC->GetPdgCode() << endl;
1545
1546                 if(TMath::Abs(negativeMC->GetPdgCode())==211 && TMath::Abs(positiveMC->GetPdgCode())==211){  // Pions ...
1547                         virtualParticleMotherLabel=virtualParticleMCLabel;
1548                         isPiPiDecay=kTRUE;
1549 //                      } else if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1550 //                              if( virtualParticleMotherMC->GetPdgCode() != 22 ){
1551 //                                      virtualParticleMotherLabel=virtualParticleMCLabel;
1552 //                                      isDalitz = kTRUE;
1553 //                              } else if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
1554 //                                      virtualParticleMotherLabel=virtualParticleMotherMC->GetFirstMother();
1555 //                                      isRealGamma = kTRUE; //no virtual gamma
1556 //                              }
1557                 }
1558         }
1559
1560         if( pi0MCLabel == virtualParticleMotherLabel ){
1561                 if(((TParticle*)MCStack->Particle(virtualParticleMotherLabel))->GetPdgCode() == 221){
1562                         isTrueEta=kTRUE;
1563                 }
1564                 if(((TParticle*)MCStack->Particle(virtualParticleMotherLabel))->GetPdgCode() == 223){
1565                         isTrueOmega=kTRUE;
1566                 }
1567         }
1568
1569         if( isTrueEta || isTrueOmega ){ // True Eta or Omega
1570                 if ( isPiPiDecay) { //real eta -> Pi+ Pi- Pi0
1571                         Float_t weighted= 1;
1572 //                      if( ((AliPrimaryPionCuts*) fPionCutArray->At(fiCut))->DoWeights() ) {
1573 //                              if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(gammaMotherLabel, fMCStack,fInputEvent)){
1574 //                                      if (((TParticle*)MCStack->Particle(gammaMotherLabel))->Pt()>0.005){
1575 //                                              weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),gammaMotherLabel,fMCStack,fInputEvent);
1576 //                                      }
1577 //                              }
1578 //                      }
1579                         fHistoTruePionPionFromNeutralMesonInvMassPt[fiCut]->Fill(TrueVirtualParticleCandidate->GetMass(),TrueVirtualParticleCandidate->Pt());
1580                         fHistoTrueMotherPiPlPiMiPiZeroInvMassPt[fiCut]->Fill(mesoncand->M(),mesoncand->Pt(),weighted);
1581                 }
1582         }
1583
1584 }
1585
1586
1587 //________________________________________________________________________
1588 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::UpdateEventByEventData(){
1589    //see header file for documentation
1590
1591         Int_t method = 1;
1592         if( method == 1 ) {
1593                 if(fNeutralPionCandidates->GetEntries() >0 ){
1594                         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
1595                                 fBGHandler[fiCut]->AddMesonEvent(fNeutralPionCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fV0Reader->GetNumberOfPrimaryTracks(),0);
1596                         } else{ // means we use #V0s for multiplicity
1597                                 fBGHandler[fiCut]->AddMesonEvent(fNeutralPionCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fGoodGammas->GetEntries(),0);
1598                         }
1599                 }
1600         } else if ( method == 2 ){
1601                 if(fGoodVirtualParticles->GetEntries() > 0 ){
1602                         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
1603                                 fBGHandler[fiCut]->AddEvent(fGoodVirtualParticles,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fV0Reader->GetNumberOfPrimaryTracks(),0);
1604                         } else{ // means we use #V0s for multiplicity
1605                                 fBGHandler[fiCut]->AddEvent(fGoodVirtualParticles,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fGoodVirtualParticles->GetEntries(),0);
1606                         }
1607                 }
1608         }
1609 }
1610
1611 //________________________________________________________________________
1612 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::MoveParticleAccordingToVertex(AliAODConversionMother* particle,const AliGammaConversionAODBGHandler::GammaConversionVertex *vertex){
1613    //see header file for documentation
1614
1615    Double_t dx = vertex->fX - fInputEvent->GetPrimaryVertex()->GetX();
1616    Double_t dy = vertex->fY - fInputEvent->GetPrimaryVertex()->GetY();
1617    Double_t dz = vertex->fZ - fInputEvent->GetPrimaryVertex()->GetZ();
1618
1619    Double_t movedPlace[3] = {particle->GetProductionX() - dx,particle->GetProductionY() - dy,particle->GetProductionZ() - dz};
1620    particle->SetProductionPoint(movedPlace);
1621 }
1622
1623 //_____________________________________________________________________________________
1624 Bool_t AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::IsEtaPiPlPiMiPiZeroDaughter( Int_t label ) const {
1625 //
1626 // Returns true if the particle comes from eta -> pi+ pi- gamma
1627 //
1628         Int_t motherLabel = fMCStack->Particle( label )->GetMother(0);
1629         if( motherLabel < 0 || motherLabel >= fMCStack->GetNtrack() ) return kFALSE;
1630
1631         TParticle* mother = fMCStack->Particle( motherLabel );
1632 //      cout << "found eta? " << endl;
1633         if( mother->GetPdgCode() != 221 ) return kFALSE;
1634 //              else cout << "YES" << endl;
1635         if( IsPiPlPiMiPiZeroDecay( mother ) ) return kTRUE;
1636         return kFALSE;
1637 }
1638
1639 //_____________________________________________________________________________________
1640 Bool_t AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::IsOmegaPiPlPiMiPiZeroDaughter( Int_t label ) const {
1641 //
1642 // Returns true if the particle comes from eta -> pi+ pi- gamma
1643 //
1644         Int_t motherLabel = fMCStack->Particle( label )->GetMother(0);
1645         if( motherLabel < 0 || motherLabel >= fMCStack->GetNtrack() ) return kFALSE;
1646
1647         TParticle* mother = fMCStack->Particle( motherLabel );
1648 //      cout << "found omega? " << endl;
1649         if( mother->GetPdgCode() != 223 ) return kFALSE;
1650 //              else cout << "YES" << endl;
1651         if( IsPiPlPiMiPiZeroDecay( mother ) ) return kTRUE;
1652         return kFALSE;
1653 }
1654
1655
1656 //_____________________________________________________________________________
1657 Bool_t AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::IsPiPlPiMiPiZeroDecay(TParticle *fMCMother) const
1658 {
1659 //      cout << fMCMother->GetNDaughters() << endl;
1660         if( fMCMother->GetNDaughters() != 3 ) return kFALSE;
1661 //      cout << fMCMother->GetPdgCode() << endl;
1662         if( !(fMCMother->GetPdgCode() == 221 || fMCMother->GetPdgCode() == 223)  ) return kFALSE;
1663 //      cout << "made it til here" << endl;
1664
1665         TParticle *posPion = 0x0;
1666         TParticle *negPion = 0x0;
1667         TParticle *neutPion    = 0x0;
1668
1669         for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
1670                 TParticle* temp = (TParticle*)fMCStack->Particle( index );
1671
1672                 switch( temp->GetPdgCode() ) {
1673                 case 211:
1674                         posPion =  temp;
1675                         break;
1676                 case -211:
1677                         negPion =  temp;
1678                         break;
1679                 case 111:
1680                         neutPion = temp;
1681                         break;
1682                 }
1683         }
1684
1685         if( posPion && negPion && neutPion) return kTRUE;
1686
1687         return kFALSE;
1688 }
1689
1690 //_____________________________________________________________________________________
1691 Bool_t AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::GammaIsNeutralMesonPiPlPiMiPiZeroDaughter( Int_t label ) const {
1692 //
1693 // Returns true if the particle comes from eta -> pi+ pi- gamma
1694 //
1695         Int_t motherLabel = fMCStack->Particle( label )->GetMother(0);
1696         if( motherLabel < 0 || motherLabel >= fMCStack->GetNtrack() ) return kFALSE;
1697
1698         TParticle* mother = fMCStack->Particle( motherLabel );
1699 //      cout << "found omega? " << endl;
1700         if( mother->GetPdgCode() != 111 ) return kFALSE;
1701 //              else cout << "YES" << endl;
1702         Int_t grandMotherLabel = mother->GetMother(0);
1703         if( grandMotherLabel < 0 || grandMotherLabel >= fMCStack->GetNtrack() ) return kFALSE;
1704         TParticle* grandmother = fMCStack->Particle( grandMotherLabel );
1705
1706         if( IsPiPlPiMiPiZeroDecay( grandmother ) ) return kTRUE;
1707         return kFALSE;
1708 }