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] = {450,250,7,4};
 261         Double_t xMin[nDim] = {0.3,0, 0,0};
 262         Double_t xMax[nDim] = {0.75,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                 }
 897         }
 898 }
 899
 900 //________________________________________________________________________
 901 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessNeutralPionCandidatesPureConversions(){
 902         // Conversion Gammas
 903         if(fGoodGammas->GetEntries()>1){
 904                 for(Int_t firstGammaIndex=0;firstGammaIndex<fGoodGammas->GetEntries()-1;firstGammaIndex++){
 905                         AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(firstGammaIndex));
 906                         if (gamma0==NULL) continue;
 907                         for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fGoodGammas->GetEntries();secondGammaIndex++){
 908                                 AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(secondGammaIndex));
 909                                 //Check for same Electron ID
 910                                 if (gamma1==NULL) continue;
 911                                 if(gamma0->GetTrackLabelPositive() == gamma1->GetTrackLabelPositive() ||
 912                                 gamma0->GetTrackLabelNegative() == gamma1->GetTrackLabelNegative() ||
 913                                 gamma0->GetTrackLabelNegative() == gamma1->GetTrackLabelPositive() ||
 914                                 gamma0->GetTrackLabelPositive() == gamma1->GetTrackLabelNegative() ) continue;
 915
 916                                 AliAODConversionMother *pi0cand = new AliAODConversionMother(gamma0,gamma1);
 917                                 pi0cand->SetLabels(firstGammaIndex,secondGammaIndex);
 918
 919                                 pi0cand->CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
 920                                 if((((AliConversionMesonCuts*)fNeutralPionMesonCutArray->At(fiCut))->MesonIsSelected(pi0cand,kTRUE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift()))){
 921                                         fHistoGammaGammaInvMassPt[fiCut]->Fill(pi0cand->M(),pi0cand->Pt());
 922                                         if (pi0cand->M() > 0.08 && pi0cand->M() < 0.145){
 923                                                 fNeutralPionCandidates->Add(pi0cand);
 924 //                                              cout << "Pi0 candidate " << pi0cand->M() << "\t" << pi0cand->Pt() << endl;
 925                                         }
 926
 927                                         if(fIsMC){
 928                                                 if(fInputEvent->IsA()==AliESDEvent::Class())
 929                                                         ProcessTrueNeutralPionCandidatesPureConversions(pi0cand,gamma0,gamma1);
 930                                                 if(fInputEvent->IsA()==AliAODEvent::Class())
 931                                                         ProcessTrueNeutralPionCandidatesPureConversionsAOD(pi0cand,gamma0,gamma1);
 932                                         }
 933                                 }
 934                                 delete pi0cand;
 935                                 pi0cand=0x0;
 936                         }
 937                 }
 938         }
 939 }
 940
 941
 942 //______________________________________________________________________
 943 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessTrueNeutralPionCandidatesPureConversions(AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1)
 944 {
 945         // Process True Mesons
 946         AliStack *MCStack = fMCEvent->Stack();
 947         if(TrueGammaCandidate0->GetV0Index()<fInputEvent->GetNumberOfV0s()){
 948                 Bool_t isTruePi0 = kFALSE;
 949                 Bool_t isTruePi0Dalitz = kFALSE;
 950                 Bool_t gamma0DalitzCand = kFALSE;
 951                 Bool_t gamma1DalitzCand = kFALSE;
 952                 Int_t gamma0MCLabel = TrueGammaCandidate0->GetMCParticleLabel(MCStack);
 953                 Int_t gamma0MotherLabel = -1;
 954                 if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
 955                         // Daughters Gamma 0
 956                         TParticle * negativeMC = (TParticle*)TrueGammaCandidate0->GetNegativeMCDaughter(MCStack);
 957                         TParticle * positiveMC = (TParticle*)TrueGammaCandidate0->GetPositiveMCDaughter(MCStack);
 958                         TParticle * gammaMC0 = (TParticle*)MCStack->Particle(gamma0MCLabel);
 959                         if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
 960                                 if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
 961                                         if(gammaMC0->GetPdgCode() == 22){ // ... with Gamma Mother
 962                                                 gamma0MotherLabel=gammaMC0->GetFirstMother();
 963                                         }
 964                                 }
 965                                 if(gammaMC0->GetPdgCode() ==111){ // Dalitz candidate
 966                                         gamma0DalitzCand = kTRUE;
 967                                         gamma0MotherLabel=-111;
 968                                 }
 969                         }
 970                 }
 971                 if(TrueGammaCandidate1->GetV0Index()<fInputEvent->GetNumberOfV0s()){
 972                         Int_t gamma1MCLabel = TrueGammaCandidate1->GetMCParticleLabel(MCStack);
 973                         Int_t gamma1MotherLabel = -1;
 974                         if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
 975                                 // Daughters Gamma 1
 976                                 TParticle * negativeMC = (TParticle*)TrueGammaCandidate1->GetNegativeMCDaughter(MCStack);
 977                                 TParticle * positiveMC = (TParticle*)TrueGammaCandidate1->GetPositiveMCDaughter(MCStack);
 978                                 TParticle * gammaMC1 = (TParticle*)MCStack->Particle(gamma1MCLabel);
 979                                 if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
 980                                         if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
 981                                                 if(gammaMC1->GetPdgCode() == 22){ // ... with Gamma Mother
 982                                                         gamma1MotherLabel=gammaMC1->GetFirstMother();
 983                                                 }
 984                                         }
 985                                         if(gammaMC1->GetPdgCode() ==111 ){ // Dalitz candidate
 986                                                 gamma1DalitzCand = kTRUE;
 987                                                 gamma1MotherLabel=-111;
 988                                         }
 989                                 }
 990                         }
 991                         if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
 992                                 if(((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetPdgCode() == 111){
 993                                         isTruePi0=kTRUE;
 994                                 }
 995                         }
 996
 997                         //Identify Dalitz candidate
 998                         if (gamma1DalitzCand || gamma0DalitzCand){
 999                                 if (gamma0DalitzCand && gamma0MCLabel >=0 && gamma0MCLabel==gamma1MotherLabel){
1000                                         if (gamma0MotherLabel == -111) isTruePi0Dalitz = kTRUE;
1001                                 }
1002                                 if (gamma1DalitzCand && gamma1MCLabel >=0 && gamma1MCLabel==gamma0MotherLabel){
1003                                         if (gamma1MotherLabel == -111) isTruePi0Dalitz = kTRUE;
1004                                 }
1005                         }
1006
1007
1008                         if(isTruePi0 || isTruePi0Dalitz){// True Pion
1009                                 fHistoTrueMotherGammaGammaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1010                         }
1011                 }
1012         }
1013 }
1014 //______________________________________________________________________
1015 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessTrueNeutralPionCandidatesPureConversionsAOD(AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1)
1016 {
1017
1018         // Process True Mesons
1019         TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1020         Bool_t isTruePi0 = kFALSE;
1021         Bool_t isTruePi0Dalitz = kFALSE;
1022         Bool_t gamma0DalitzCand = kFALSE;
1023         Bool_t gamma1DalitzCand = kFALSE;
1024
1025         if (AODMCTrackArray!=NULL && TrueGammaCandidate0 != NULL){
1026                 AliAODMCParticle *positiveMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate0->GetMCLabelPositive()));
1027                 AliAODMCParticle *negativeMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate0->GetMCLabelNegative()));
1028
1029                 Int_t gamma0MCLabel = -1;
1030                 Int_t gamma0MotherLabel = -1;
1031                 if(!positiveMC||!negativeMC)
1032                         return;
1033
1034                 if(positiveMC->GetMother()>-1&&(negativeMC->GetMother() == positiveMC->GetMother())){
1035                         gamma0MCLabel = positiveMC->GetMother();
1036                 }
1037
1038                 if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1039                         // Daughters Gamma 0
1040                         AliAODMCParticle * gammaMC0 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma0MCLabel));
1041                         if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1042                                 if(((positiveMC->GetMCProcessCode())) == 5 && ((negativeMC->GetMCProcessCode())) == 5){ // ... From Conversion ...
1043                                         if(gammaMC0->GetPdgCode() == 22){ // ... with Gamma Mother
1044                                         gamma0MotherLabel=gammaMC0->GetMother();
1045                                         }
1046                                 }
1047                                 if(gammaMC0->GetPdgCode() ==111){ // Dalitz candidate
1048                                         gamma0DalitzCand = kTRUE;
1049                                         gamma0MotherLabel=-111;
1050                                 }
1051                         }
1052                 }
1053                 positiveMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate1->GetMCLabelPositive()));
1054                 negativeMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate1->GetMCLabelNegative()));
1055
1056                 Int_t gamma1MCLabel = -1;
1057                 Int_t gamma1MotherLabel = -1;
1058                 if(!positiveMC||!negativeMC)
1059                         return;
1060
1061                 if(positiveMC->GetMother()>-1&&(negativeMC->GetMother() == positiveMC->GetMother())){
1062                         gamma1MCLabel = positiveMC->GetMother();
1063                 }
1064                 if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1065                         // Daughters Gamma 1
1066                         AliAODMCParticle * gammaMC1 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MCLabel));
1067                         if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1068                                 if(((positiveMC->GetMCProcessCode())) == 5 && ((negativeMC->GetMCProcessCode())) == 5){ // ... From Conversion ...
1069                                         if(gammaMC1->GetPdgCode() == 22){ // ... with Gamma Mother
1070                                         gamma1MotherLabel=gammaMC1->GetMother();
1071                                         }
1072                                 }
1073                                 if(gammaMC1->GetPdgCode() ==111 ){ // Dalitz candidate
1074                                                 gamma1DalitzCand = kTRUE;
1075                                                 gamma1MotherLabel=-111;
1076                                 }
1077                         }
1078                 }
1079                 if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
1080                         if(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->GetPdgCode() == 111){
1081                                 isTruePi0=kTRUE;
1082                         }
1083                 }
1084
1085                 //Identify Dalitz candidate
1086                 if (gamma1DalitzCand || gamma0DalitzCand){
1087                         if (gamma0DalitzCand && gamma0MCLabel >=0 && gamma0MCLabel==gamma1MotherLabel){
1088                                 if (gamma0MotherLabel == -111) isTruePi0Dalitz = kTRUE;
1089                         }
1090                         if (gamma1DalitzCand && gamma1MCLabel >=0 && gamma1MCLabel==gamma0MotherLabel){
1091                                 if (gamma1MotherLabel == -111) isTruePi0Dalitz = kTRUE;
1092                         }
1093                 }
1094
1095                 if(isTruePi0 || isTruePi0Dalitz){// True Pion
1096                         fHistoTrueMotherGammaGammaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1097                 }
1098         }
1099         return;
1100 }
1101
1102
1103
1104 //________________________________________________________________________
1105 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessPionCandidates(){
1106
1107         Double_t magField = fInputEvent->GetMagneticField();
1108         if( magField  < 0.0 ){
1109                 magField =  1.0;
1110         } else {
1111                 magField =  -1.0;
1112         }
1113
1114         vector<Int_t> lGoodNegPionIndexPrev(0);
1115         vector<Int_t> lGoodPosPionIndexPrev(0);
1116
1117         for(UInt_t i = 0; i < fSelectorNegPionIndex.size(); i++){
1118                 AliESDtrack* negPionCandidate = fESDEvent->GetTrack(fSelectorNegPionIndex[i]);
1119                 if(! ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->PionIsSelected(negPionCandidate) ) continue;
1120                 lGoodNegPionIndexPrev.push_back(   fSelectorNegPionIndex[i] );
1121                 fHistoNegPionPt[fiCut]->Fill(negPionCandidate->Pt());
1122                 fHistoNegPionPhi[fiCut]->Fill(negPionCandidate->Phi());
1123                 if( fMCEvent ) {
1124                         Int_t labelNegPion = TMath::Abs( negPionCandidate->GetLabel() );
1125                         if( labelNegPion < fMCStack->GetNtrack() ){
1126                                 TParticle* negPion = fMCStack->Particle(labelNegPion);
1127                                 if( negPion->GetPdgCode() ==  -211 ){
1128                                         if( labelNegPion < fMCStack->GetNprimary() ){
1129                                                 fHistoTrueNegPionPt[fiCut]->Fill(negPionCandidate->Pt());    //primary negPion
1130                                         }
1131                                         if( IsEtaPiPlPiMiPiZeroDaughter(labelNegPion) || IsOmegaPiPlPiMiPiZeroDaughter(labelNegPion) ) {
1132                                                 if( labelNegPion < fMCStack->GetNprimary() ) {
1133                                                         fHistoTrueNegPionFromNeutralMesonPt[fiCut]->Fill(negPionCandidate->Pt());
1134                                                 }
1135                                         }
1136                                 }
1137                         }
1138                 }
1139         }
1140
1141         for(UInt_t i = 0; i < fSelectorPosPionIndex.size(); i++){
1142                 AliESDtrack* posPionCandidate = fESDEvent->GetTrack( fSelectorPosPionIndex[i] );
1143                 if(! ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->PionIsSelected(posPionCandidate) ) continue;
1144                 lGoodPosPionIndexPrev.push_back(   fSelectorPosPionIndex[i]  );
1145                 fHistoPosPionPt[fiCut]->Fill( posPionCandidate->Pt() );
1146                 fHistoPosPionPhi[fiCut]->Fill( posPionCandidate->Phi() );
1147
1148                 if( fMCEvent ) {
1149                         Int_t labelPosPion = TMath::Abs( posPionCandidate->GetLabel() );
1150                         if( labelPosPion < fMCStack->GetNtrack() ) {
1151                                 TParticle* posPion = fMCStack->Particle(labelPosPion);
1152                                 if( posPion->GetPdgCode() ==  211 ){
1153                                         if( labelPosPion < fMCStack->GetNprimary() ){
1154                                                 fHistoTruePosPionPt[fiCut]->Fill(posPionCandidate->Pt());
1155                                         }
1156                                         if( IsEtaPiPlPiMiPiZeroDaughter(labelPosPion) || IsOmegaPiPlPiMiPiZeroDaughter(labelPosPion) ) {
1157                                                 if( labelPosPion < fMCStack->GetNprimary() ){
1158                                                         fHistoTruePosPionFromNeutralMesonPt[fiCut]->Fill(posPionCandidate->Pt());
1159                                                 }
1160                                         }
1161                                 }
1162                         }
1163                 }
1164         }
1165
1166
1167         for(UInt_t i = 0; i < lGoodNegPionIndexPrev.size(); i++){
1168
1169                 AliESDtrack *negPionCandidate = fESDEvent->GetTrack(lGoodNegPionIndexPrev[i]);
1170                 AliKFParticle negPionCandidateKF( *negPionCandidate->GetConstrainedParam(), 211 );
1171
1172                 for(UInt_t j = 0; j < lGoodPosPionIndexPrev.size(); j++){
1173
1174                         AliESDtrack *posPionCandidate = fESDEvent->GetTrack(lGoodPosPionIndexPrev[j]);
1175                         AliKFParticle posPionCandidateKF( *posPionCandidate->GetConstrainedParam(), 211 );
1176
1177                         AliKFConversionPhoton* virtualPhoton = NULL;
1178                         virtualPhoton = new AliKFConversionPhoton(negPionCandidateKF,posPionCandidateKF);
1179                         AliKFVertex primaryVertexImproved(*fInputEvent->GetPrimaryVertex());
1180 //                      primaryVertexImproved+=*virtualPhoton;
1181                         virtualPhoton->SetProductionVertex(primaryVertexImproved);
1182                         virtualPhoton->SetTrackLabels( lGoodPosPionIndexPrev[j], lGoodNegPionIndexPrev[i]);
1183
1184
1185                         if( fMCEvent ) {
1186                                 Int_t labeln=TMath::Abs(negPionCandidate->GetLabel());
1187                                 Int_t labelp=TMath::Abs(posPionCandidate->GetLabel());
1188                                 TParticle *fNegativeMCParticle = fMCStack->Particle(labeln);
1189                                 TParticle *fPositiveMCParticle = fMCStack->Particle(labelp);
1190
1191                                 if( fPositiveMCParticle && fNegativeMCParticle) {
1192                                         virtualPhoton->SetMCLabelPositive(labelp);
1193                                         virtualPhoton->SetMCLabelNegative(labeln);
1194                                 }
1195                         }
1196
1197                         AliAODConversionPhoton *vParticle = new AliAODConversionPhoton(virtualPhoton); //To Apply PsiPairCut
1198                         if (((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->DoMassCut()){
1199                                 if (vParticle->GetMass() < ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->GetMassCut()){
1200                                         fGoodVirtualParticles->Add(  vParticle );
1201                                 }
1202                         } else {
1203                                 fGoodVirtualParticles->Add(  vParticle );
1204                         }
1205                         delete virtualPhoton;
1206                         virtualPhoton=NULL;
1207
1208                 }
1209         }
1210 }
1211
1212 //_____________________________________________________________________________
1213 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessMCParticles(){
1214
1215         // Loop over all primary MC particle
1216
1217         for(Int_t i = 0; i < fMCStack->GetNprimary(); i++) {
1218
1219                 TParticle* particle = (TParticle *)fMCStack->Particle(i);
1220                 if (!particle) continue;
1221
1222                 Int_t isMCFromMBHeader = -1;
1223                 if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 0){
1224                         isMCFromMBHeader
1225                                 = ((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent);
1226                         if(isMCFromMBHeader == 0 && ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1227                 }
1228
1229                 if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack,fInputEvent)){
1230
1231                         // find MC photons
1232                         if(((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->PhotonIsSelectedMC(particle,fMCStack,kFALSE)){
1233                                 fHistoMCAllGammaPt[fiCut]->Fill(particle->Pt()); // All MC Gamma
1234                                 if(particle->GetMother(0) >-1){
1235                                         if (fMCStack->Particle(particle->GetMother(0))->GetPdgCode() ==111){
1236                                                 if (fMCStack->Particle(particle->GetMother(0))->GetMother(0) > -1){
1237                                                         if ( fMCStack->Particle((fMCStack->Particle(particle->GetMother(0)))->GetMother(0))->GetPdgCode() == 221 ||
1238                                                                  fMCStack->Particle((fMCStack->Particle(particle->GetMother(0)))->GetMother(0))->GetPdgCode() == 223 ){
1239                                                                 if ( fMCStack->Particle(particle->GetMother(0))->GetNDaughters()==3 )
1240                                                                         fHistoMCGammaFromNeutralMesonPt[fiCut]->Fill(particle->Pt()); // All photons from eta or omega via pi0
1241                                                         }
1242                                                 }
1243                                         }
1244                                 }
1245                         }
1246
1247                         if(((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->PhotonIsSelectedMC(particle,fMCStack,kTRUE)){
1248                                 fHistoMCConvGammaPt[fiCut]->Fill(particle->Pt());
1249                         } // Converted MC Gamma
1250
1251                         if(((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->PionIsSelectedMC(i,fMCStack)){
1252                                 if( particle->GetPdgCode() == 211){
1253                                         fHistoMCAllPosPionsPt[fiCut]->Fill(particle->Pt()); // All pos pions
1254                                         if(particle->GetMother(0) >-1){
1255                                                 if (fMCStack->Particle(particle->GetMother(0))->GetPdgCode() ==221 || fMCStack->Particle(particle->GetMother(0))->GetPdgCode() ==223)
1256                                                         fHistoMCPosPionsFromNeutralMesonPt[fiCut]->Fill(particle->Pt()); // All pos from eta or omega
1257                                         }
1258                                 }
1259                                 if( particle->GetPdgCode() == -211){
1260                                         fHistoMCAllNegPionsPt[fiCut]->Fill(particle->Pt()); // All neg pions
1261                                         if(particle->GetMother(0) >-1){
1262                                                 if (fMCStack->Particle(particle->GetMother(0))->GetPdgCode() ==221 || fMCStack->Particle(particle->GetMother(0))->GetPdgCode() ==223 )
1263                                                         fHistoMCNegPionsFromNeutralMesonPt[fiCut]->Fill(particle->Pt()); // All pos from eta or omega
1264                                         }
1265                                 }
1266                         }
1267
1268
1269                         // \eta -> pi+ pi- \gamma
1270                         Int_t labelNeutPion = -1;
1271                         Int_t labelNegPion = -1;
1272                         Int_t labelPosPion = -1;
1273
1274                         if( ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelectedMCPiPlPiMiPiZero(particle,fMCStack,labelNegPion,labelPosPion,labelNeutPion,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift())){
1275                                 Float_t weighted= 1;
1276                                 if( ((AliPrimaryPionCuts*) fPionCutArray->At(fiCut))->DoWeights() ) {
1277                                         if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack,fInputEvent)){
1278                                                 if (particle->Pt()>0.005){
1279                                                         weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, fMCStack,fInputEvent);
1280                                                 }
1281                                         }
1282                                 }
1283                                 if(particle->GetPdgCode() == 221)fHistoMCEtaPiPlPiMiPiZeroPt[fiCut]->Fill(particle->Pt(), weighted);                                            // All MC Eta in respective decay channel
1284                                 if(particle->GetPdgCode() == 223)fHistoMCOmegaPiPlPiMiPiZeroPt[fiCut]->Fill(particle->Pt(), weighted);                                          // All MC Omega in respective decay channel
1285
1286                                 TParticle *neutPion    = fMCStack->Particle(labelNeutPion);
1287                                 TParticle *gamma1 = fMCStack->Particle(neutPion->GetDaughter(0));
1288                                 TParticle *gamma2 = fMCStack->Particle(neutPion->GetDaughter(1));
1289                                 if(
1290                                         ((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->PhotonIsSelectedMC(gamma1,fMCStack,kFALSE) &&                                    // test first daugther of pi0
1291                                         ((AliConversionPhotonCuts*)fGammaCutArray->At(fiCut))->PhotonIsSelectedMC(gamma2,fMCStack,kFALSE) &&                                    // test second daughter of pi0
1292                                         ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->PionIsSelectedMC(labelNegPion,fMCStack) &&                                                             // test negative pion
1293                                         ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->PionIsSelectedMC(labelPosPion,fMCStack)                                                                // test positive pion
1294                                    ) {
1295                                                 if(particle->GetPdgCode() == 221) fHistoMCEtaPiPlPiMiPiZeroInAccPt[fiCut]->Fill(particle->Pt(), weighted );             // MC Eta pi+ pi- pi0 with gamma's and e+e- in acc
1296                                                 if(particle->GetPdgCode() == 223) fHistoMCOmegaPiPlPiMiPiZeroInAccPt[fiCut]->Fill(particle->Pt(), weighted );           // MC Omega pi+ pi- pi0 with gamma's and e+e- in acc
1297                                 }
1298                         }
1299                 }
1300         }
1301 }
1302
1303
1304 //________________________________________________________________________
1305 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::CalculateMesonCandidates(){
1306
1307         // Conversion Gammas
1308         if( fNeutralPionCandidates->GetEntries() > 0 && fGoodVirtualParticles->GetEntries() > 0 ){
1309
1310                 vector<Bool_t> lGoodVirtualParticle(fGoodVirtualParticles->GetEntries(), kFALSE);
1311
1312                 for(Int_t mesonIndex=0; mesonIndex<fNeutralPionCandidates->GetEntries(); mesonIndex++){
1313                         AliAODConversionMother *neutralPion=static_cast<AliAODConversionMother*>(fNeutralPionCandidates->At(mesonIndex));
1314                         if (neutralPion==NULL) continue;
1315                         for(Int_t virtualParticleIndex=0;virtualParticleIndex<fGoodVirtualParticles->GetEntries();virtualParticleIndex++){
1316
1317                                 AliAODConversionPhoton *vParticle=dynamic_cast<AliAODConversionPhoton*>(fGoodVirtualParticles->At(virtualParticleIndex));
1318                                 if (vParticle==NULL) continue;
1319                                 //Check for same Electron ID
1320
1321                                 AliAODConversionMother *mesoncand = new AliAODConversionMother(neutralPion,vParticle);
1322                                 mesoncand->SetLabels(mesonIndex,virtualParticleIndex);
1323
1324                                 if( ( ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(mesoncand,kTRUE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift())) ){
1325
1326 //                                      cout<< "Meson Accepted "<<endl;
1327
1328                                         Int_t zbin= fBGHandler[fiCut]->GetZBinIndex(fESDEvent->GetPrimaryVertex()->GetZ());
1329                                         Int_t mbin = 0;
1330                                         if( ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
1331                                                 mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fNumberOfESDTracks);
1332                                         } else {
1333                                                 mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGoodGammas->GetEntries());
1334                                         }
1335
1336                                         AliESDtrack *posPionVParticle = 0;
1337                                         AliESDtrack *negPionVParticle = 0;
1338
1339                                         Double_t clsToFPos = -1.0;
1340                                         Double_t clsToFNeg = -1.0;
1341
1342                                         Float_t dcaToVertexXYPos = -1.0;
1343                                         Float_t dcaToVertexZPos  = -1.0;
1344                                         Float_t dcaToVertexXYNeg = -1.0;
1345                                         Float_t dcaToVertexZNeg  = -1.0;
1346
1347
1348                                         if ( fDoMesonQA ) {
1349
1350                                                 fHistoPionPionInvMassPt[fiCut]->Fill( vParticle->GetMass(),vParticle->Pt());
1351
1352                                                 posPionVParticle = fESDEvent->GetTrack( vParticle->GetTrackLabelPositive() );
1353                                                 negPionVParticle = fESDEvent->GetTrack( vParticle->GetTrackLabelNegative() );
1354                                                 clsToFPos = ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->GetNFindableClustersTPC(posPionVParticle);
1355                                                 clsToFNeg = ((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->GetNFindableClustersTPC(negPionVParticle);
1356
1357                                                 Float_t bPos[2];
1358                                                 Float_t bCovPos[3];
1359                                                 posPionVParticle->GetImpactParameters(bPos,bCovPos);
1360                                                 if (bCovPos[0]<=0 || bCovPos[2]<=0) {
1361                                                         AliDebug(1, "Estimated b resolution lower or equal zero!");
1362                                                         bCovPos[0]=0; bCovPos[2]=0;
1363                                                 }
1364
1365                                                 Float_t bNeg[2];
1366                                                 Float_t bCovNeg[3];
1367                                                 posPionVParticle->GetImpactParameters(bNeg,bCovNeg);
1368                                                 if (bCovNeg[0]<=0 || bCovNeg[2]<=0) {
1369                                                         AliDebug(1, "Estimated b resolution lower or equal zero!");
1370                                                         bCovNeg[0]=0; bCovNeg[2]=0;
1371                                                 }
1372
1373                                                 dcaToVertexXYPos = bPos[0];
1374                                                 dcaToVertexZPos  = bPos[1];
1375                                                 dcaToVertexXYNeg = bNeg[0];
1376                                                 dcaToVertexZNeg  = bNeg[1];
1377                                         }
1378
1379                                         fHistoMotherInvMassPt[fiCut]->Fill(mesoncand->M(),mesoncand->Pt());
1380                                         Double_t sparesFill[4] = {mesoncand->M(),mesoncand->Pt(),(Double_t)zbin,(Double_t)mbin};
1381                                         fTHnSparseMotherInvMassPtZM[fiCut]->Fill(sparesFill,1);
1382
1383                                         if ( fDoMesonQA ) {
1384                                                 if( lGoodVirtualParticle[virtualParticleIndex] == kFALSE ) {
1385
1386                                                         fHistoNegPionEta[fiCut]->Fill( negPionVParticle->Eta() );
1387                                                         fHistoPosPionEta[fiCut]->Fill( posPionVParticle->Eta() );
1388
1389                                                         fHistoNegPionClsTPC[fiCut]->Fill(clsToFNeg,negPionVParticle->Pt());
1390                                                         fHistoPosPionClsTPC[fiCut]->Fill(clsToFPos,posPionVParticle->Pt());
1391
1392                                                         fHistoPionDCAxy[fiCut]->Fill(  dcaToVertexXYNeg, negPionVParticle->Pt() );
1393                                                         fHistoPionDCAz[fiCut]->Fill(   dcaToVertexZNeg,  negPionVParticle->Pt() );
1394                                                         fHistoPionDCAxy[fiCut]->Fill(  dcaToVertexXYPos, posPionVParticle->Pt() );
1395                                                         fHistoPionDCAz[fiCut]->Fill(   dcaToVertexZPos,  posPionVParticle->Pt() );
1396
1397                                                         fHistoPionTPCdEdxNSigma[fiCut]->Fill( posPionVParticle->P(),((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->GetPIDResponse()->NumberOfSigmasTPC(posPionVParticle, AliPID::kPion) );
1398                                                         fHistoPionTPCdEdxNSigma[fiCut]->Fill( negPionVParticle->P(),((AliPrimaryPionCuts*)fPionCutArray->At(fiCut))->GetPIDResponse()->NumberOfSigmasTPC(negPionVParticle, AliPID::kPion) );
1399
1400                                                         fHistoPionTPCdEdx[fiCut]->Fill( posPionVParticle->P(), TMath::Abs(posPionVParticle->GetTPCsignal()));
1401                                                         fHistoPionTPCdEdx[fiCut]->Fill( negPionVParticle->P(), TMath::Abs(negPionVParticle->GetTPCsignal()));
1402
1403                                                         lGoodVirtualParticle[virtualParticleIndex] = kTRUE;
1404
1405                                                 }
1406                                         }
1407
1408
1409                                         if(fMCEvent){
1410                                                 ProcessTrueMesonCandidates(mesoncand,neutralPion,vParticle);
1411                                         }
1412                                 }
1413                                 delete mesoncand;
1414                                 mesoncand=0x0;
1415                         }
1416                 }
1417         }
1418 }
1419
1420 //________________________________________________________________________
1421 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::CalculateBackground(){
1422
1423         Int_t zbin= fBGHandler[fiCut]->GetZBinIndex(fESDEvent->GetPrimaryVertex()->GetZ());
1424         Int_t mbin = 0;
1425
1426
1427         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
1428                 mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fNumberOfESDTracks);
1429         } else {
1430                 mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGoodGammas->GetEntries());
1431         }
1432
1433         Int_t method = 1;
1434         AliGammaConversionAODBGHandler::GammaConversionVertex *bgEventVertex = NULL;
1435         if( ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity() ) {
1436                 for(Int_t nEventsInBG=0;nEventsInBG<fBGHandler[fiCut]->GetNBGEvents();nEventsInBG++){
1437                         AliGammaConversionAODVector *previousEventV0s = fBGHandler[fiCut]->GetBGGoodV0s(zbin,mbin,nEventsInBG);
1438                         if(fMoveParticleAccordingToVertex == kTRUE && method == 1){
1439                                 bgEventVertex = fBGHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG);
1440                         }
1441
1442                         for(Int_t iCurrent=0;iCurrent<fGoodVirtualParticles->GetEntries();iCurrent++){
1443                                 AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGoodVirtualParticles->At(iCurrent));
1444
1445                                 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
1446                                         AliAODConversionPhoton previousGoodV0 = (AliAODConversionPhoton)(*(previousEventV0s->at(iPrevious)));
1447
1448                                         if(fMoveParticleAccordingToVertex == kTRUE && method == 1 ){
1449                                                 MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
1450                                         }
1451
1452                                         AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&previousGoodV0);
1453
1454
1455                                         if( ( ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE, ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift()))){
1456                                                 fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
1457                                                 Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin};
1458                                                 fTHnSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
1459                                         }
1460                                         delete backgroundCandidate;
1461                                         backgroundCandidate = 0x0;
1462                                 }
1463                         }
1464                 }
1465         } else {
1466                 for(Int_t nEventsInBG=0;nEventsInBG <fBGHandler[fiCut]->GetNBGEvents();nEventsInBG++){
1467                         AliGammaConversionAODVector *previousEventV0s = fBGHandler[fiCut]->GetBGGoodV0s(zbin,mbin,nEventsInBG);
1468                         if(previousEventV0s){
1469                                 if(fMoveParticleAccordingToVertex == kTRUE && method == 1){
1470                                         bgEventVertex = fBGHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG);
1471                                 }
1472                                 for(Int_t iCurrent=0;iCurrent<fGoodVirtualParticles->GetEntries();iCurrent++){
1473
1474                                         AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGoodVirtualParticles->At(iCurrent));
1475
1476                                         for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
1477
1478                                                 AliAODConversionPhoton previousGoodV0 = (AliAODConversionPhoton)(*(previousEventV0s->at(iPrevious)));
1479
1480                                                 if(fMoveParticleAccordingToVertex == kTRUE && method ==1){
1481                                                         MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
1482                                                 }
1483
1484                                                 AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&previousGoodV0);
1485
1486                                                 if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE,((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetEtaShift()))){
1487                                                         fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
1488                                                         Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin};
1489                                                         fTHnSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
1490                                                 }
1491                                                 delete backgroundCandidate;
1492                                                 backgroundCandidate = 0x0;
1493                                         }
1494                                 }
1495                         }
1496                 }
1497         }
1498 }
1499
1500 //______________________________________________________________________
1501 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::ProcessTrueMesonCandidates(AliAODConversionMother *mesoncand, AliAODConversionMother *TrueNeutralPionCandidate, AliAODConversionPhoton *TrueVirtualParticleCandidate){
1502
1503 //      // Process True Mesons
1504 //
1505 //      AliStack *MCStack = fMCEvent->Stack();
1506 //
1507 //      Bool_t isTrueEta = kFALSE;
1508 //      Bool_t isTrueOmega = kFALSE;
1509 //      Int_t gammaMCLabel = TrueNeutralPionCandidate->GetMCParticleLabel(MCStack);
1510 //      Int_t gammaMotherLabel = -1;
1511 //
1512 //      if(gammaMCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1513 //              // Daughters Gamma 0
1514 //              TParticle * negativeMC = (TParticle*)TrueNeutralPionCandidate->GetNegativeMCDaughter(MCStack);
1515 //              TParticle * positiveMC = (TParticle*)TrueNeutralPionCandidate->GetPositiveMCDaughter(MCStack);
1516 //              TParticle * gammaMC = (TParticle*)MCStack->Particle(gammaMCLabel);
1517 //
1518 //              if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1519 //                      if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
1520 //                              if(gammaMC->GetPdgCode() == 22){ // ... with Gamma Mother
1521 //                                      gammaMotherLabel=gammaMC->GetFirstMother();
1522 //                              }
1523 //                      }
1524 //              }
1525 //
1526 //
1527 //      }
1528 //
1529 //      Int_t virtualParticleMCLabel = TrueVirtualParticleCandidate->GetMCParticleLabel(MCStack);
1530 //      Int_t virtualParticleMotherLabel = -1;
1531 //
1532 //      Bool_t isPiPiDecay = kFALSE;
1533 //
1534 //      if (fDoMesonQA){
1535 //              TParticle * negativeMC = (TParticle*)TrueVirtualParticleCandidate->GetNegativeMCDaughter(MCStack);
1536 //              TParticle * positiveMC = (TParticle*)TrueVirtualParticleCandidate->GetPositiveMCDaughter(MCStack);
1537 //              if(TMath::Abs(negativeMC->GetPdgCode())==211 && TMath::Abs(positiveMC->GetPdgCode())==211){  // Pions ...
1538 //                      fHistoTruePionPionInvMassPt[fiCut]->Fill(TrueVirtualParticleCandidate->GetMass(),TrueVirtualParticleCandidate->Pt());
1539 //              }
1540 //      }
1541 //
1542 //      if(virtualParticleMCLabel != -1){ // if virtualParticleMCLabel==-1 particles don't have same mother
1543 //              TParticle * negativeMC = (TParticle*)TrueVirtualParticleCandidate->GetNegativeMCDaughter(MCStack);
1544 //              TParticle * positiveMC = (TParticle*)TrueVirtualParticleCandidate->GetPositiveMCDaughter(MCStack);
1545 // //                   TParticle * virtualParticleMotherMC = (TParticle*)MCStack->Particle(virtualParticleMCLabel);
1546 // //                   cout << "pdg code same mother - " << virtualParticleMotherMC->GetPdgCode() << endl;
1547 //
1548 //              if(TMath::Abs(negativeMC->GetPdgCode())==211 && TMath::Abs(positiveMC->GetPdgCode())==211){  // Pions ...
1549 //                      virtualParticleMotherLabel=virtualParticleMCLabel;
1550 //                      isPiPiDecay=kTRUE;
1551 // //                   } else if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1552 // //                           if( virtualParticleMotherMC->GetPdgCode() != 22 ){
1553 // //                                   virtualParticleMotherLabel=virtualParticleMCLabel;
1554 // //                                   isDalitz = kTRUE;
1555 // //                           } else if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
1556 // //                                   virtualParticleMotherLabel=virtualParticleMotherMC->GetFirstMother();
1557 // //                                   isRealGamma = kTRUE; //no virtual gamma
1558 // //                           }
1559 //              }
1560 //      }
1561 //
1562 //      if(gammaMotherLabel >= 0 && ( gammaMotherLabel == virtualParticleMotherLabel) ){
1563 //              if(((TParticle*)MCStack->Particle(virtualParticleMotherLabel))->GetPdgCode() == 221){
1564 //                      isTrueEta=kTRUE;
1565 //              }
1566 //              if(((TParticle*)MCStack->Particle(virtualParticleMotherLabel))->GetPdgCode() == 223){
1567 //                      isTrueOmega=kTRUE;
1568 //              }
1569 //      }
1570 //
1571 //      if( isTrueEta || isTrueOmega ){ // True Eta or Omega
1572 //              if ( isPiPiDecay) { //real eta -> Pi+ Pi- Pi0
1573 //                      Float_t weighted= 1;
1574 //                      if( ((AliPrimaryPionCuts*) fPionCutArray->At(fiCut))->DoWeights() ) {
1575 //                              if(((AliConvEventCuts*)fEventCutArray->At(fiCut))->IsParticleFromBGEvent(gammaMotherLabel, fMCStack,fInputEvent)){
1576 //                                      if (((TParticle*)MCStack->Particle(gammaMotherLabel))->Pt()>0.005){
1577 //                                              weighted= ((AliConvEventCuts*)fEventCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),gammaMotherLabel,fMCStack,fInputEvent);
1578 //                                      }
1579 //                              }
1580 //                      }
1581 //                      fHistoTruePionPionFromNeutralMesonInvMassPt[fiCut]->Fill(TrueVirtualParticleCandidate->GetMass(),TrueVirtualParticleCandidate->Pt());
1582 //                      fHistoTrueMotherPiPlPiMiPiZeroInvMassPt[fiCut]->Fill(EtaCandidate->M(),EtaCandidate->Pt(),weighted);
1583 //              }
1584 //      }
1585
1586 }
1587
1588
1589 //________________________________________________________________________
1590 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::UpdateEventByEventData(){
1591    //see header file for documentation
1592
1593         Int_t method = 1;
1594         if( method == 1 ) {
1595                 if(fGoodGammas->GetEntries() >0 ){
1596                         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
1597                                 fBGHandler[fiCut]->AddEvent(fGoodGammas,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fV0Reader->GetNumberOfPrimaryTracks(),0);
1598                         } else{ // means we use #V0s for multiplicity
1599                                 fBGHandler[fiCut]->AddEvent(fGoodGammas,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fGoodGammas->GetEntries(),0);
1600                         }
1601                 }
1602         } else if ( method == 2 ){
1603                 if(fGoodVirtualParticles->GetEntries() > 0 ){
1604                         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
1605                                 fBGHandler[fiCut]->AddEvent(fGoodVirtualParticles,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fV0Reader->GetNumberOfPrimaryTracks(),0);
1606                         } else{ // means we use #V0s for multiplicity
1607                                 fBGHandler[fiCut]->AddEvent(fGoodVirtualParticles,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fGoodVirtualParticles->GetEntries(),0);
1608                         }
1609                 }
1610         }
1611 }
1612
1613 //________________________________________________________________________
1614 void AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::MoveParticleAccordingToVertex(AliAODConversionPhoton* particle,const AliGammaConversionAODBGHandler::GammaConversionVertex *vertex){
1615    //see header file for documentation
1616
1617    Double_t dx = vertex->fX - fInputEvent->GetPrimaryVertex()->GetX();
1618    Double_t dy = vertex->fY - fInputEvent->GetPrimaryVertex()->GetY();
1619    Double_t dz = vertex->fZ - fInputEvent->GetPrimaryVertex()->GetZ();
1620
1621    Double_t movedPlace[3] = {particle->GetConversionX() - dx,particle->GetConversionY() - dy,particle->GetConversionZ() - dz};
1622    particle->SetConversionPoint(movedPlace);
1623 }
1624
1625 //_____________________________________________________________________________________
1626 Bool_t AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::IsEtaPiPlPiMiPiZeroDaughter( Int_t label ) const {
1627 //
1628 // Returns true if the particle comes from eta -> pi+ pi- gamma
1629 //
1630         Int_t motherLabel = fMCStack->Particle( label )->GetMother(0);
1631         if( motherLabel < 0 || motherLabel >= fMCStack->GetNtrack() ) return kFALSE;
1632
1633         TParticle* mother = fMCStack->Particle( motherLabel );
1634         if( mother->GetPdgCode() != 221 ) return kFALSE;
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         if( mother->GetPdgCode() != 223 ) return kFALSE;
1649         if( IsPiPlPiMiPiZeroDecay( mother ) ) return kTRUE;
1650         return kFALSE;
1651 }
1652
1653
1654 //_____________________________________________________________________________
1655 Bool_t AliAnalysisTaskNeutralMesonToPiPlPiMiPiZero::IsPiPlPiMiPiZeroDecay(TParticle *fMCMother) const
1656 {
1657
1658         if( fMCMother->GetNDaughters() != 3 ) return kFALSE;
1659         if( fMCMother->GetPdgCode() != 221 || fMCMother->GetPdgCode() != 223  ) return kFALSE;
1660
1661         TParticle *posPion = 0x0;
1662         TParticle *negPion = 0x0;
1663         TParticle *neutPion    = 0x0;
1664
1665         for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
1666                 TParticle* temp = (TParticle*)fMCStack->Particle( index );
1667
1668                 switch( temp->GetPdgCode() ) {
1669                 case 211:
1670                         posPion =  temp;
1671                         break;
1672                 case -211:
1673                         negPion =  temp;
1674                         break;
1675                 case 111:
1676                         neutPion = temp;
1677                         break;
1678                 }
1679         }
1680
1681         if( posPion && negPion && neutPion) return kTRUE;
1682
1683         return kFALSE;
1684 }