PWGGA/GammaConv/AliAnalysisTaskGammaConvDalitzV1.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: Pedro González, Pedro Ladrón de Guevara, Ernesto López Torres, *
   5  *         Eulogio Serradilla                                             *
   6  * Version 2                                                           *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 // Analysis task for pi0->e+e-gamma (Dalitz decay)
  18
  19 #include <vector>
  20
  21 #include "TParticle.h"
  22 #include "TPDGCode.h"
  23 #include "TMCProcess.h"
  24 #include "TDatabasePDG.h"
  25 #include "TList.h"
  26 #include "TChain.h"
  27 #include "TDirectory.h"
  28 #include "TTree.h"
  29 #include "TH1.h"
  30 #include "TH1F.h"
  31 #include "THnSparse.h"
  32 #include "TH2F.h"
  33 #include "AliStack.h"
  34 #include "AliAnalysisManager.h"
  35 #include "AliESDInputHandler.h"
  36 #include "AliESDtrack.h"
  37 #include "AliMCEvent.h"
  38 #include "AliStack.h"
  39 #include "AliMCEventHandler.h"
  40 #include "AliPID.h"
  41 #include "AliLog.h"
  42 #include "AliESDtrackCuts.h"
  43 #include "AliESDpidCuts.h"
  44 #include "AliMCEvent.h"
  45 #include "AliESDv0.h"
  46 #include "AliESDEvent.h"
  47 #include "AliESDpid.h"
  48 #include "AliKFParticle.h"
  49 #include "AliMCEventHandler.h"
  50 #include "AliKFVertex.h"
  51 #include "AliTriggerAnalysis.h"
  52 #include "AliCentrality.h"
  53 #include "AliMultiplicity.h"
  54 #include "AliAnalysisTaskGammaConvDalitzV1.h"
  55
  56
  57 ClassImp( AliAnalysisTaskGammaConvDalitzV1 )
  58
  59 //-----------------------------------------------------------------------------------------------
  60 AliAnalysisTaskGammaConvDalitzV1::AliAnalysisTaskGammaConvDalitzV1():
  61 fV0Reader(NULL),
  62    fElecSelector(NULL),
  63    fBGHandler(NULL),
  64    fESDEvent(NULL),
  65    fMCEvent(NULL),
  66    fMCStack(NULL),
  67    fCutFolder(NULL),
  68    fESDList(NULL),
  69    fBackList(NULL),
  70    fTrueList(NULL),
  71    fMCList(NULL),
  72    fOutputContainer(0),
  73    fReaderGammas(NULL),
  74    fSelectorElectronIndex(0),
  75    fSelectorPositronIndex(0),
  76    fGoodGammas(NULL),
  77    fGoodVirtualGammas(NULL),
  78    fGoodElectrons(NULL),
  79    fGoodPositrons(NULL),
  80    fCutGammaArray(NULL),
  81    fCutElectronArray(NULL),
  82    fCutMesonArray(NULL),
  83    fGammasPool(NULL),
  84    fConversionCuts(NULL),
  85    hESDConvGammaPt(NULL),
  86    hESDDalitzElectronPt(NULL),
  87    hESDDalitzPositronPt(NULL),
  88    hESDEposEnegPsiPairDPhi(NULL),
  89    hESDEposEnegInvMassPt(NULL),
  90    hESDEposEnegLikeSignBackInvMassPt(NULL),
  91    hESDMotherInvMassPt(NULL),
  92    hESDPi0MotherInvMassPt(NULL),
  93    hESDPi0MotherDiffInvMassPt(NULL),
  94    sESDMotherInvMassPtZM(NULL),
  95    hESDMotherBackInvMassPt(NULL),
  96    sESDMotherBackInvMassPtZM(NULL),
  97    hMCAllGammaPt(NULL),
  98    hMCAllPositronsPt(NULL),
  99    hMCAllElectronsPt(NULL),
 100    hMCPi0DalitzGammaPt(NULL),
 101    hMCPi0DalitzElectronPt(NULL),
 102    hMCPi0DalitzPositronPt(NULL),
 103    hMCPi0Pt(NULL),
 104    hMCPi0GGPt(NULL),
 105    hMCEtaPt(NULL),
 106    hMCEtaGGPt(NULL),
 107    hMCPi0InAccPt(NULL),
 108    hMCEtaInAccPt(NULL),
 109    hESDEposEnegTruePi0DalitzInvMassPt(NULL),
 110    hESDEposEnegTrueEtaDalitzInvMassPt(NULL),
 111    hESDEposEnegTruePhotonInvMassPt(NULL),
 112    hESDTrueMotherInvMassPt(NULL),
 113    hESDTrueMotherPi0GGInvMassPt(NULL),
 114    hESDTruePrimaryMotherInvMassMCPt(NULL),
 115    hESDTruePrimaryPi0DalitzESDPtMCPt(NULL),
 116    hESDTrueSecondaryMotherInvMassPt(NULL),
 117    hESDTrueSecondaryMotherFromK0sInvMassPt(NULL),
 118    hESDTrueBckGGInvMassPt(NULL),
 119    hESDTrueBckContInvMassPt(NULL),
 120    hESDTrueMotherGGInvMassPt(NULL),
 121    hESDTrueConvGammaPt(NULL),
 122    hESDTruePositronPt(NULL),
 123    hESDTrueElectronPt(NULL),
 124    hESDTruePi0DalitzConvGammaPt(NULL),
 125    hESDTruePi0DalitzPositronPt(NULL),
 126    hESDTruePi0DalitzElectronPt(NULL),
 127    hNEvents(NULL),
 128    hNGoodESDTracks(NULL),
 129    fRandom(0),
 130    fUnsmearedPx(NULL),
 131    fUnsmearedPy(NULL),
 132    fUnsmearedPz(NULL),
 133    fUnsmearedE(NULL),
 134    fnCuts(0),
 135    fiCut(0),
 136    fNumberOfESDTracks(0),
 137    fMoveParticleAccordingToVertex(kFALSE),
 138    fIsHeavyIon(kFALSE),
 139    fDoMesonAnalysis(kTRUE)
 140 {
 141
 142 }
 143
 144 //-----------------------------------------------------------------------------------------------
 145 AliAnalysisTaskGammaConvDalitzV1::AliAnalysisTaskGammaConvDalitzV1( const char* name ):
 146    AliAnalysisTaskSE(name),
 147    fV0Reader(NULL),
 148    fElecSelector(NULL),
 149    fBGHandler(NULL),
 150    fESDEvent(NULL),
 151    fMCEvent(NULL),
 152    fMCStack(NULL),
 153    fCutFolder(NULL),
 154    fESDList(NULL),
 155    fBackList(NULL),
 156    fTrueList(NULL),
 157    fMCList(NULL),
 158    fOutputContainer(0),
 159    fReaderGammas(NULL),
 160    fSelectorElectronIndex(0),
 161    fSelectorPositronIndex(0),
 162    fGoodGammas(NULL),
 163    fGoodVirtualGammas(NULL),
 164    fGoodElectrons(NULL),
 165    fGoodPositrons(NULL),
 166    fCutGammaArray(NULL),
 167    fCutElectronArray(NULL),
 168    fCutMesonArray(NULL),
 169    fGammasPool(NULL),
 170    fConversionCuts(NULL),
 171    hESDConvGammaPt(NULL),
 172    hESDDalitzElectronPt(NULL),
 173    hESDDalitzPositronPt(NULL),
 174    hESDEposEnegPsiPairDPhi(NULL),
 175    hESDEposEnegInvMassPt(NULL),
 176    hESDEposEnegLikeSignBackInvMassPt(NULL),
 177    hESDMotherInvMassPt(NULL),
 178    hESDPi0MotherInvMassPt(NULL),
 179    hESDPi0MotherDiffInvMassPt(NULL),
 180    sESDMotherInvMassPtZM(NULL),
 181    hESDMotherBackInvMassPt(NULL),
 182    sESDMotherBackInvMassPtZM(NULL),
 183    hMCAllGammaPt(NULL),
 184    hMCAllPositronsPt(NULL),
 185    hMCAllElectronsPt(NULL),
 186    hMCPi0DalitzGammaPt(NULL),
 187    hMCPi0DalitzElectronPt(NULL),
 188    hMCPi0DalitzPositronPt(NULL),
 189    hMCPi0Pt(NULL),
 190    hMCPi0GGPt(NULL),
 191    hMCEtaPt(NULL),
 192    hMCEtaGGPt(NULL),
 193    hMCPi0InAccPt(NULL),
 194    hMCEtaInAccPt(NULL),
 195    hESDEposEnegTruePi0DalitzInvMassPt(NULL),
 196    hESDEposEnegTrueEtaDalitzInvMassPt(NULL),
 197    hESDEposEnegTruePhotonInvMassPt(NULL),
 198    hESDTrueMotherInvMassPt(NULL),
 199    hESDTrueMotherPi0GGInvMassPt(NULL),
 200    hESDTruePrimaryMotherInvMassMCPt(NULL),
 201    hESDTruePrimaryPi0DalitzESDPtMCPt(NULL),
 202    hESDTrueSecondaryMotherInvMassPt(NULL),
 203    hESDTrueSecondaryMotherFromK0sInvMassPt(NULL),
 204    hESDTrueBckGGInvMassPt(NULL),
 205    hESDTrueBckContInvMassPt(NULL),
 206    hESDTrueMotherGGInvMassPt(NULL),
 207    hESDTrueConvGammaPt(NULL),
 208    hESDTruePositronPt(NULL),
 209    hESDTrueElectronPt(NULL),
 210    hESDTruePi0DalitzConvGammaPt(NULL),
 211    hESDTruePi0DalitzPositronPt(NULL),
 212    hESDTruePi0DalitzElectronPt(NULL),
 213    hNEvents(NULL),
 214    hNGoodESDTracks(NULL),
 215    fRandom(0),
 216    fUnsmearedPx(NULL),
 217    fUnsmearedPy(NULL),
 218    fUnsmearedPz(NULL),
 219    fUnsmearedE(NULL),
 220    fnCuts(0),
 221    fiCut(0),
 222    fNumberOfESDTracks(0),
 223    fMoveParticleAccordingToVertex(kFALSE),
 224    fIsHeavyIon(kFALSE),
 225    fDoMesonAnalysis(kTRUE)
 226 {
 227    DefineOutput(1, TList::Class());
 228 }
 229
 230 //-----------------------------------------------------------------------------------------------
 231 AliAnalysisTaskGammaConvDalitzV1::~AliAnalysisTaskGammaConvDalitzV1()
 232 {
 233    //
 234    // virtual destructor
 235    //
 236    cout<<"Destructor"<<endl;
 237
 238    if(fGoodGammas){
 239       delete fGoodGammas;
 240       fGoodGammas = 0x0;
 241    }
 242    if(fGoodVirtualGammas){
 243       delete fGoodVirtualGammas;
 244       fGoodGammas = 0x0;
 245    }
 246    if(fGoodElectrons){
 247       delete fGoodGammas;
 248       fGoodGammas = 0x0;
 249    }
 250    if(fGoodPositrons){
 251       delete fGoodGammas;
 252       fGoodGammas = 0x0;
 253    }
 254    if(fBGHandler){
 255       delete[] fBGHandler;
 256       fBGHandler = 0x0;
 257    }
 258    if( fGammasPool ){
 259       delete[] fGammasPool;
 260       fGammasPool = 0x0;
 261    }
 262 }
 263
 264 //___________________________________________________________
 265 void AliAnalysisTaskGammaConvDalitzV1::InitBack(){
 266
 267    Double_t *zBinLimitsArray= new Double_t[9];
 268    zBinLimitsArray[0] = -50.00;
 269    zBinLimitsArray[1] = -3.375;
 270    zBinLimitsArray[2] = -1.605;
 271    zBinLimitsArray[3] = -0.225;
 272    zBinLimitsArray[4] = 1.065;
 273    zBinLimitsArray[5] = 2.445;
 274    zBinLimitsArray[6] = 4.245;
 275    zBinLimitsArray[7] = 50.00;
 276    zBinLimitsArray[8] = 1000.00;
 277
 278    Double_t *multiplicityBinLimitsArrayTracks= new Double_t[6];
 279    multiplicityBinLimitsArrayTracks[0] = 0;
 280    multiplicityBinLimitsArrayTracks[1] = 8.5;
 281    multiplicityBinLimitsArrayTracks[2] = 16.5;
 282    multiplicityBinLimitsArrayTracks[3] = 27.5;
 283    multiplicityBinLimitsArrayTracks[4] = 41.5;
 284    multiplicityBinLimitsArrayTracks[5] = 200.;
 285
 286    if(fIsHeavyIon){
 287       multiplicityBinLimitsArrayTracks[0] = 0;
 288       multiplicityBinLimitsArrayTracks[1] = 200.;
 289       multiplicityBinLimitsArrayTracks[2] = 500.;
 290       multiplicityBinLimitsArrayTracks[3] = 1000.;
 291       multiplicityBinLimitsArrayTracks[4] = 1500.;
 292       multiplicityBinLimitsArrayTracks[5] = 5000.;
 293    }
 294
 295    Double_t *multiplicityBinLimitsArrayV0s= new Double_t[5];
 296
 297    multiplicityBinLimitsArrayV0s[0] = 2;
 298    multiplicityBinLimitsArrayV0s[1] = 3;
 299    multiplicityBinLimitsArrayV0s[2] = 4;
 300    multiplicityBinLimitsArrayV0s[3] = 5;
 301    multiplicityBinLimitsArrayV0s[4] = 9999;
 302
 303    if(fIsHeavyIon){
 304       multiplicityBinLimitsArrayV0s[0] = 2;
 305       multiplicityBinLimitsArrayV0s[1] = 10;
 306       multiplicityBinLimitsArrayV0s[2] = 30;
 307       multiplicityBinLimitsArrayV0s[3] = 50;
 308       multiplicityBinLimitsArrayV0s[4] = 9999;
 309    }
 310
 311    const Int_t nDim = 4;
 312    Int_t    nBins[nDim] = {1000,250,8,5};
 313    Double_t xMin[nDim]  = {0,0, 0,0};
 314    Double_t xMax[nDim]  = {1,25,8,5};
 315
 316    sESDMotherInvMassPtZM     = new THnSparseF*[fnCuts];
 317    sESDMotherBackInvMassPtZM = new THnSparseF*[fnCuts];
 318
 319
 320    fBGHandler = new AliGammaConversionAODBGHandler*[fnCuts];
 321
 322    for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 323
 324
 325       TString cutstringElectron     =   ((AliDalitzElectronCuts*)fCutElectronArray->At(iCut))->GetCutNumber();
 326       TString cutstringMeson        =   ((AliConversionMesonCuts*)fCutMesonArray->At(iCut))->GetCutNumber();
 327       TString cutstringGamma        =   ((AliConversionCuts*)fCutGammaArray->At(iCut))->GetCutNumber();
 328
 329
 330       fBackList[iCut] = new TList();
 331       fBackList[iCut]->SetName(Form("%s_%s_%s Back histograms",cutstringGamma.Data(),cutstringElectron.Data(),cutstringMeson.Data()));
 332       fBackList[iCut]->SetOwner(kTRUE);
 333       fCutFolder[iCut]->Add(fBackList[iCut]);
 334
 335       sESDMotherInvMassPtZM[iCut] = new THnSparseF("Back_Mother_InvMass_Pt_z_m","Back_Mother_InvMass_Pt_z_m",nDim,nBins,xMin,xMax);
 336       sESDMotherInvMassPtZM[iCut]->Sumw2();
 337       fBackList[iCut]->Add(sESDMotherInvMassPtZM[iCut]);
 338       sESDMotherBackInvMassPtZM[iCut] = new THnSparseF("Back_Back_InvMass_Pt_z_m","Back_Back_InvMass_Pt_z_m",nDim,nBins,xMin,xMax);
 339       sESDMotherBackInvMassPtZM[iCut]->Sumw2();
 340       fBackList[iCut]->Add(sESDMotherBackInvMassPtZM[iCut]);
 341
 342       if(((AliDalitzElectronCuts*)fCutElectronArray->At(iCut))->UseTrackMultiplicity()) {
 343          fBGHandler[iCut] = new AliGammaConversionAODBGHandler(9,6,((AliDalitzElectronCuts*)fCutElectronArray->At(iCut))->NumberOfRotationEvents());
 344          fBGHandler[iCut]->Initialize(zBinLimitsArray, multiplicityBinLimitsArrayTracks);
 345       }
 346       else{
 347          fBGHandler[iCut] = new AliGammaConversionAODBGHandler(9,5,((AliDalitzElectronCuts*)fCutElectronArray->At(iCut))->NumberOfRotationEvents());
 348          fBGHandler[iCut]->Initialize(zBinLimitsArray, multiplicityBinLimitsArrayV0s);
 349       }
 350       if( ( (AliDalitzElectronCuts*)fCutElectronArray->At(iCut))->GetBKGMethod() == 3 ){
 351          fGammasPool[iCut] = new TList();
 352       }
 353    }
 354
 355 }
 356
 357 //______________________________________________________________________
 358 void AliAnalysisTaskGammaConvDalitzV1::UserCreateOutputObjects()
 359 {
 360    //
 361    // Create ouput objects
 362    //
 363
 364    // Create the output container
 365    if(fOutputContainer != NULL){
 366       delete fOutputContainer;
 367       fOutputContainer = NULL;
 368    }
 369    if(fOutputContainer == NULL){
 370       fOutputContainer = new TList();
 371       fOutputContainer->SetOwner(kTRUE);
 372    }
 373
 374    fGoodGammas = new TList();
 375    //fGoodGammas->SetOwner(kTRUE);
 376
 377
 378    fGoodVirtualGammas = new TList();
 379    //fGoodVirtualGammas->SetOwner(kTRUE);
 380
 381
 382
 383    fGammasPool                     = new TList*[fnCuts];
 384    fCutFolder                      = new TList*[fnCuts];
 385    fESDList                        = new TList*[fnCuts];
 386    fBackList                       = new TList*[fnCuts];
 387    hNEvents                        = new TH1I*[fnCuts];
 388    hNGoodESDTracks                 = new TH1I*[fnCuts];
 389    hESDConvGammaPt                 = new TH1F*[fnCuts];
 390    hESDDalitzElectronPt            = new TH1F*[fnCuts];
 391    hESDDalitzPositronPt            = new TH1F*[fnCuts];
 392    hESDEposEnegPsiPairDPhi         = new TH2F*[fnCuts];
 393    hESDEposEnegInvMassPt           = new TH2F*[fnCuts];
 394    hESDEposEnegLikeSignBackInvMassPt = new TH2F*[fnCuts];
 395    hESDMotherInvMassPt             = new TH2F*[fnCuts];
 396    hESDPi0MotherInvMassPt          = new TH2F*[fnCuts];
 397    hESDPi0MotherDiffInvMassPt      = new TH2F*[fnCuts];
 398    hESDMotherBackInvMassPt         = new TH2F*[fnCuts];
 399
 400
 401    for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 402
 403
 404       TString cutstringElectron =((AliDalitzElectronCuts*)fCutElectronArray->At(iCut))->GetCutNumber();
 405       TString cutstringMeson= ((AliConversionMesonCuts*)fCutMesonArray->At(iCut))->GetCutNumber();
 406       TString cutstringGamma = ((AliConversionCuts*)fCutGammaArray->At(iCut))->GetCutNumber();
 407
 408       fCutFolder[iCut] = new TList();
 409       fCutFolder[iCut]->SetName(Form("Cut Number %s_%s_%s",cutstringGamma.Data(),cutstringElectron.Data(),cutstringMeson.Data()));
 410       fCutFolder[iCut]->SetOwner(kTRUE);
 411       fOutputContainer->Add(fCutFolder[iCut]);
 412
 413       fESDList[iCut] = new TList();
 414       fESDList[iCut]->SetName(Form("%s_%s_%s ESD histograms",cutstringGamma.Data(),cutstringElectron.Data(),cutstringMeson.Data()));
 415       fESDList[iCut]->SetOwner(kTRUE);
 416
 417
 418
 419       hNEvents[iCut] = new TH1I("NEvents","NEvents",9,-0.5,8.5);
 420       hNEvents[iCut]->GetXaxis()->SetBinLabel(1,"Accepted");
 421       hNEvents[iCut]->GetXaxis()->SetBinLabel(2,"Centrality");
 422       hNEvents[iCut]->GetXaxis()->SetBinLabel(3,"Missing MC");
 423       hNEvents[iCut]->GetXaxis()->SetBinLabel(4,"Trigger");
 424       hNEvents[iCut]->GetXaxis()->SetBinLabel(5,"Vertex Z");
 425       hNEvents[iCut]->GetXaxis()->SetBinLabel(6,"Cont. Vertex");
 426       hNEvents[iCut]->GetXaxis()->SetBinLabel(7,"Pile-Up");
 427       hNEvents[iCut]->GetXaxis()->SetBinLabel(8,"no SDD");
 428       hNEvents[iCut]->GetXaxis()->SetBinLabel(9,"no V0AND");
 429       fESDList[iCut]->Add(hNEvents[iCut]);
 430
 431
 432
 433       if(fIsHeavyIon) hNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",3000,0,3000);
 434       else hNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",200,0,200);
 435       fESDList[iCut]->Add(hNGoodESDTracks[iCut]);
 436
 437       hESDConvGammaPt[iCut] = new TH1F("ESD_ConvGamma_Pt","ESD_ConvGamma_Pt",250,0,25);
 438       fESDList[iCut]->Add(hESDConvGammaPt[iCut]);
 439
 440       hESDDalitzElectronPt[iCut] = new TH1F("ESD_DalitzElectron_Pt","ESD_DalitzElectron_Pt",250,0,25);
 441       fESDList[iCut]->Add(hESDDalitzElectronPt[iCut]);
 442
 443       hESDDalitzPositronPt[iCut] = new TH1F("ESD_DalitzPositron_Pt","ESD_DalitzPositron_Pt",250,0,25);
 444       fESDList[iCut]->Add(hESDDalitzPositronPt[iCut]);
 445
 446
 447       hESDEposEnegPsiPairDPhi[iCut] = new TH2F("ESD_EposEneg_PsiPair_DPhi","ESD_EposEneg_PsiPair_DPhi", 100, -1.0,1.0,100,-1.0,1.0 );
 448       fESDList[iCut]->Add(hESDEposEnegPsiPairDPhi[iCut]);
 449
 450       hESDEposEnegInvMassPt[iCut] = new TH2F("ESD_EposEneg_InvMassPt","ESD_EposEneg_InvMassPt",5000,0.,5.,100,0.,10.);
 451       fESDList[iCut]->Add(hESDEposEnegInvMassPt[iCut]);
 452
 453       hESDEposEnegLikeSignBackInvMassPt[iCut]  = new TH2F("ESD_EposEneg_LikeSignBack_InvMassPt","ESD_EposEneg_LikeSignBack_InvMassPt",5000,0.,5.,100,0.,10.);
 454       fESDList[iCut]->Add(hESDEposEnegLikeSignBackInvMassPt[iCut]);
 455
 456
 457       hESDMotherInvMassPt[iCut] = new TH2F("ESD_DalitzMother_InvMass_Pt","ESD_DalitzMother_InvMass_Pt",1000,0,1,250,0,25);
 458       fESDList[iCut]->Add(hESDMotherInvMassPt[iCut]);
 459
 460       hESDPi0MotherInvMassPt[iCut] = new TH2F("ESD_Pi0Mother_InvMass_Pt","ESD_Pi0Mother_InvMass_Pt",4000,0,4,250,0,25);
 461       fESDList[iCut]->Add(hESDPi0MotherInvMassPt[iCut]);
 462
 463       hESDPi0MotherDiffInvMassPt[iCut] = new TH2F("ESD_Pi0Mother_DiffInvMass_Pt","ESD_Pi0Mother_DiffInvMass_Pt",2000,0,2,250,0,25);
 464       fESDList[iCut]->Add(hESDPi0MotherDiffInvMassPt[iCut]);
 465
 466
 467       hESDMotherBackInvMassPt[iCut] = new TH2F("ESD_DalitzBackground_InvMass_Pt","ESD_DalitzBackground_InvMass_Pt",1000,0,1,250,0,25);
 468       fESDList[iCut]->Add(hESDMotherBackInvMassPt[iCut]);
 469
 470       fCutFolder[iCut]->Add(fESDList[iCut]);
 471
 472
 473
 474    }
 475
 476
 477    InitBack(); // Init Background Handler
 478
 479
 480    if(MCEvent()){
 481       // MC Histogramms
 482       fMCList = new TList*[fnCuts];
 483       // True Histogramms
 484       fTrueList = new TList*[fnCuts];
 485       hESDTrueConvGammaPt = new TH1F*[fnCuts];
 486       hESDTruePositronPt  = new TH1F*[fnCuts];
 487       hESDTrueElectronPt  = new TH1F*[fnCuts];
 488       hESDTruePi0DalitzConvGammaPt = new TH1F*[fnCuts];
 489       hESDTruePi0DalitzPositronPt  = new TH1F*[fnCuts];
 490       hESDTruePi0DalitzElectronPt  = new TH1F*[fnCuts];
 491       //if(fDoMesonAnalysis){
 492       hMCAllGammaPt = new TH1F*[fnCuts];
 493       hMCAllPositronsPt = new TH1F*[fnCuts];
 494       hMCAllElectronsPt = new TH1F*[fnCuts];
 495       hMCPi0DalitzGammaPt    = new TH1F*[fnCuts];
 496       hMCPi0DalitzElectronPt = new TH1F*[fnCuts];
 497       hMCPi0DalitzPositronPt = new TH1F*[fnCuts];
 498
 499       hMCPi0Pt = new TH1F*[fnCuts];
 500       hMCPi0GGPt =  new TH1F*[fnCuts];
 501       hMCEtaPt = new TH1F*[fnCuts];
 502       hMCEtaGGPt = new TH1F*[fnCuts];
 503       hMCPi0InAccPt = new TH1F*[fnCuts];
 504       hMCEtaInAccPt = new TH1F*[fnCuts];
 505
 506
 507       hESDEposEnegTruePi0DalitzInvMassPt           = new TH2F*[fnCuts];
 508       hESDEposEnegTrueEtaDalitzInvMassPt           = new TH2F*[fnCuts];
 509       hESDEposEnegTruePhotonInvMassPt              = new TH2F*[fnCuts];
 510
 511       hESDTrueMotherInvMassPt = new TH2F*[fnCuts];
 512       hESDTrueMotherPi0GGInvMassPt = new TH2F*[fnCuts];
 513       hESDTruePrimaryPi0DalitzESDPtMCPt = new TH2F*[fnCuts];
 514       hESDTruePrimaryMotherInvMassMCPt = new TH2F*[fnCuts];
 515       hESDTrueSecondaryMotherInvMassPt = new TH2F*[fnCuts];
 516       hESDTrueSecondaryMotherFromK0sInvMassPt = new TH2F*[fnCuts];
 517       hESDTrueBckGGInvMassPt = new TH2F*[fnCuts];
 518       hESDTrueBckContInvMassPt = new TH2F*[fnCuts];
 519       hESDTrueMotherGGInvMassPt = new TH2F*[fnCuts];
 520       //}
 521
 522       for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 523          TString cutstringElectron =((AliDalitzElectronCuts*)fCutElectronArray->At(iCut))->GetCutNumber();
 524          TString cutstringMeson= ((AliConversionMesonCuts*)fCutMesonArray->At(iCut))->GetCutNumber();
 525          TString cutstringGamma = ((AliConversionCuts*)fCutGammaArray->At(iCut))->GetCutNumber();
 526
 527          fMCList[iCut] = new TList();
 528          fMCList[iCut]->SetName(Form("%s_%s_%s MC histograms",cutstringGamma.Data(),cutstringElectron.Data(),cutstringMeson.Data()));
 529          fMCList[iCut]->SetOwner(kTRUE);
 530          fCutFolder[iCut]->Add(fMCList[iCut]);
 531
 532
 533                 hMCAllGammaPt[iCut] = new TH1F("MC_AllGamma_Pt","MC_AllGamma_Pt",250,0,25);
 534                 fMCList[iCut]->Add(hMCAllGammaPt[iCut]);
 535
 536                 hMCAllPositronsPt[iCut] = new TH1F("MC_AllPositrons_Pt","MC_AllPositrons_Pt",250,0,25);
 537                 fMCList[iCut]->Add(hMCAllPositronsPt[iCut]);
 538
 539                 hMCAllElectronsPt[iCut] = new TH1F("MC_AllElectrons_Pt","MC_AllElectrons_Pt",250,0,25);
 540                 fMCList[iCut]->Add(hMCAllElectronsPt[iCut]);
 541
 542                 hMCPi0DalitzGammaPt[iCut] = new TH1F("MC_Pi0DalitzGamma_Pt","MC_Pi0DalitzGamma_Pt",250,0,25);
 543                 fMCList[iCut]->Add(hMCPi0DalitzGammaPt[iCut]);
 544
 545                 hMCPi0DalitzPositronPt[iCut] = new TH1F("MC_Pi0DalitzPositron_Pt","MC_Pi0DalitzPositron_Pt",250,0,25);
 546                 fMCList[iCut]->Add(hMCPi0DalitzPositronPt[iCut]);
 547
 548                 hMCPi0DalitzElectronPt[iCut] = new TH1F("MC_Pi0DalitzElectron_Pt","MC_Pi0DalitzElectron_Pt",250,0,25);
 549                 fMCList[iCut]->Add(hMCPi0DalitzElectronPt[iCut]);
 550
 551
 552                 hMCPi0Pt[iCut] = new TH1F("MC_Pi0_Pt","MC_Pi0_Pt",250,0,25);
 553                 fMCList[iCut]->Add(hMCPi0Pt[iCut]);
 554
 555                 hMCPi0GGPt[iCut] = new TH1F("MC_Pi0_GG_Pt","MC_Pi0_GG_Pt",250,0,25);
 556                 fMCList[iCut]->Add(hMCPi0GGPt[iCut]);
 557
 558
 559                 hMCEtaPt[iCut] = new TH1F("MC_Eta_Pt","MC_Eta_Pt",250,0,25);
 560                 fMCList[iCut]->Add(hMCEtaPt[iCut]);
 561
 562                 hMCEtaGGPt[iCut] = new TH1F("MC_Eta_GG_Pt","MC_Eta_GG_Pt",250,0,25);
 563                 fMCList[iCut]->Add(hMCEtaGGPt[iCut]);
 564
 565
 566                 hMCPi0InAccPt[iCut] = new TH1F("MC_Pi0DalitzInAcc_Pt","MC_Pi0DalitzInAcc_Pt",250,0,25);
 567                 fMCList[iCut]->Add(hMCPi0InAccPt[iCut]);
 568                 hMCEtaInAccPt[iCut] = new TH1F("MC_EtaDalitzInAcc_Pt","MC_EtaDalitzInAcc_Pt",250,0,25);
 569                 fMCList[iCut]->Add(hMCEtaInAccPt[iCut]);
 570
 571          fTrueList[iCut] = new TList();
 572          fTrueList[iCut]->SetName(Form("%s_%s_%s True histograms",cutstringGamma.Data(),cutstringElectron.Data(),cutstringMeson.Data()));
 573          fTrueList[iCut]->SetOwner(kTRUE);
 574          fCutFolder[iCut]->Add(fTrueList[iCut]);
 575
 576
 577
 578
 579
 580
 581
 582
 583
 584
 585          hESDEposEnegTruePi0DalitzInvMassPt[iCut] = new TH2F("ESD_EposEneg_TruePi0Dalitz_InvMassPt","ESD_EposEneg_TruePi0Dalitz_InvMassPt",5000,0.,5.,100,0.,10.);
 586          fTrueList[iCut]->Add(hESDEposEnegTruePi0DalitzInvMassPt[iCut]);
 587
 588          hESDEposEnegTrueEtaDalitzInvMassPt[iCut] = new TH2F("ESD_EposEneg_TrueEtaDalitz_InvMassPt","ESD_EposEneg_TrueEtaDalitz_InvMassPt",5000,0.,5.,100,0.,10.);
 589          fTrueList[iCut]->Add(hESDEposEnegTrueEtaDalitzInvMassPt[iCut]);
 590
 591          hESDEposEnegTruePhotonInvMassPt[iCut] = new TH2F("ESD_EposEneg_TruePhoton_InvMassPt","ESD_EposEneg_TruePhoton_InvMassPt",5000,0.,5.,100,0.,10.);
 592          fTrueList[iCut]->Add(hESDEposEnegTruePhotonInvMassPt[iCut]);
 593
 594
 595
 596          hESDTruePositronPt[iCut] = new TH1F("ESD_TruePositron_Pt","ESD_TruePositron_Pt",250,0,25);
 597          fTrueList[iCut]->Add(hESDTruePositronPt[iCut]);
 598
 599          hESDTrueElectronPt[iCut] = new TH1F("ESD_TrueElectron_Pt","ESD_TrueElectron_Pt",250,0,25);
 600          fTrueList[iCut]->Add(hESDTrueElectronPt[iCut]);
 601
 602          hESDTrueConvGammaPt[iCut] = new TH1F("ESD_TrueConvGamma_Pt","ESD_TrueConvGamma_Pt",250,0,25);
 603          fTrueList[iCut]->Add(hESDTrueConvGammaPt[iCut]);
 604
 605          hESDTruePi0DalitzConvGammaPt[iCut] = new TH1F("ESD_TruePi0DalitzConvGamma_Pt","ESD_TruePi0DalitzConvGamma_Pt",250,0,25);
 606          fTrueList[iCut]->Add(hESDTruePi0DalitzConvGammaPt[iCut]);
 607
 608          hESDTruePi0DalitzElectronPt[iCut] = new TH1F("ESD_TruePi0DalitzElectron_Pt","ESD_TruePi0DalitzElectron_Pt",250,0,25);
 609          fTrueList[iCut]->Add(hESDTruePi0DalitzElectronPt[iCut]);
 610
 611          hESDTruePi0DalitzPositronPt[iCut] = new TH1F("ESD_TruePi0DalitzPositron_Pt","ESD_TruePi0DalitzPositron_Pt",250,0,25);
 612          fTrueList[iCut]->Add(hESDTruePi0DalitzPositronPt[iCut]);
 613
 614
 615
 616
 617
 618          hESDTrueMotherInvMassPt[iCut] = new TH2F("ESD_TrueMother_InvMass_Pt","ESD_TrueMother_InvMass_Pt",1000,0,1,250,0,25);
 619          fTrueList[iCut]->Add(hESDTrueMotherInvMassPt[iCut]);
 620
 621          hESDTrueMotherPi0GGInvMassPt[iCut] = new TH2F("ESD_TrueMotherPi0GG_InvMass_Pt","ESD_TrueMotherPi0GG_InvMass_Pt",1000,0,1,250,0,25);
 622          fTrueList[iCut]->Add(hESDTrueMotherPi0GGInvMassPt[iCut]);
 623          hESDTruePrimaryPi0DalitzESDPtMCPt[iCut] = new TH2F("ESD_TruePrimaryPi0Dalitz_ESDPt_MCPt","ESD_TruePrimaryPi0Dalitz_ESDPt_MCPt",250,0,25,250,0,25);
 624          fTrueList[iCut]->Add(hESDTruePrimaryPi0DalitzESDPtMCPt[iCut]);
 625          hESDTruePrimaryMotherInvMassMCPt[iCut] = new TH2F("ESD_TruePrimaryMother_InvMass_MCPt","ESD_TrueDalitzPrimaryMother_InvMass_MCPt",1000,0,1,250,0,25);
 626          fTrueList[iCut]->Add(hESDTruePrimaryMotherInvMassMCPt[iCut]);
 627          hESDTrueSecondaryMotherInvMassPt[iCut] = new TH2F("ESD_TrueDalitzSecondaryMother_InvMass_Pt","ESD_TrueDalitzSecondaryMother_InvMass_Pt",1000,0,1,250,0,25);
 628          fTrueList[iCut]->Add(hESDTrueSecondaryMotherInvMassPt[iCut]);
 629          //                             hESDTrueSecondaryMotherFromK0sInvMassPt[iCut] = new TH2F("ESD_TrueDalitzSecondaryMotherFromK0s_InvMass_Pt","ESD_TrueDalitzSecondaryMotherFromK0s_InvMass_Pt",1000,0,1,250,0,25);
 630          //                             fTrueList[iCut]->Add(hESDTrueSecondaryMotherFromK0sInvMassPt[iCut]);
 631          hESDTrueBckGGInvMassPt[iCut] = new TH2F("ESD_TrueDalitzBckGG_InvMass_Pt","ESD_TrueDalitzBckGG_InvMass_Pt",1000,0,1,250,0,25);
 632          fTrueList[iCut]->Add(hESDTrueBckGGInvMassPt[iCut]);
 633          hESDTrueBckContInvMassPt[iCut] = new TH2F("ESD_TrueDalitzBckCont_InvMass_Pt","ESD_TrueDalitzBckCont_InvMass_Pt",1000,0,1,250,0,25);
 634          fTrueList[iCut]->Add(hESDTrueBckContInvMassPt[iCut]);
 635          //                             hESDTrueMotherGGInvMassPt[iCut] = new TH2F("ESD_TrueGammaGamma_InvMass_Pt","ESD_TrueGammaGamma_InvMass_Pt",1000,0,1,250,0,25);
 636          //                             fTrueList[iCut]->Add(hESDTrueMotherGGInvMassPt[iCut]);
 637
 638       }
 639    }
 640
 641    PostData(1, fOutputContainer);
 642
 643 }
 644
 645 //______________________________________________________________________
 646 void AliAnalysisTaskGammaConvDalitzV1::UserExec(Option_t *)
 647 {
 648
 649    //
 650    // Execute analysis for current event
 651    //
 652
 653
 654    fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
 655    if(!fV0Reader){printf("Error: No V0 Reader");return;} // GetV0Reader
 656
 657
 658    Int_t eventQuality = ((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetEventQuality();
 659
 660    if(eventQuality == 2 || eventQuality == 3){// Event Not Accepted due to MC event missing or wrong trigger for V0ReaderV1
 661       for(Int_t iCut = 0; iCut<fnCuts; iCut++){
 662          hNEvents[iCut]->Fill(eventQuality);
 663       }
 664       return;
 665    }
 666
 667
 668
 669
 670    fElecSelector=(AliDalitzElectronSelector*)AliAnalysisManager::GetAnalysisManager()->GetTask("ElectronSelector");
 671    if(!fElecSelector){printf("Error: No ElectronSelector");return;} // GetV0Reader
 672
 673
 674    fMCEvent         =  MCEvent();
 675    fESDEvent        = (AliESDEvent*)InputEvent();
 676    fReaderGammas    = fV0Reader->GetReconstructedGammas(); // Gammas from default Cut
 677    fSelectorElectronIndex = fElecSelector->GetReconstructedElectronsIndex(); // Electrons from default Cut
 678    fSelectorPositronIndex = fElecSelector->GetReconstructedPositronsIndex(); // Positrons from default Cut
 679
 680    CountESDTracks(); // Estimate Event Multiplicity
 681    //AddTaskContainers(); //Add conatiner
 682
 683    for(Int_t iCut = 0; iCut<fnCuts; iCut++){
 684       fiCut = iCut;
 685
 686       Int_t eventNotAccepted =
 687          ((AliConversionCuts*)fCutGammaArray->At(iCut))
 688          ->IsEventAcceptedByConversionCut(fV0Reader->GetConversionCuts(),fInputEvent,fMCEvent,fIsHeavyIon);
 689       if(eventNotAccepted){
 690          //                     cout << "event rejected due to wrong trigger: " <<eventNotAccepted << endl;
 691          hNEvents[iCut]->Fill(eventNotAccepted); // Check Centrality, PileUp, SDD and V0AND --> Not Accepted => eventQuality = 1
 692          continue;
 693       }
 694
 695       if(eventQuality != 0){// Event Not Accepted
 696          //                     cout << "event rejected due to: " <<eventQuality << endl;
 697          hNEvents[iCut]->Fill(eventQuality);
 698          continue;
 699       }
 700
 701       hNEvents[iCut]->Fill(eventQuality);
 702
 703       hNGoodESDTracks[iCut]->Fill(fNumberOfESDTracks);
 704
 705       if(fMCEvent){ // Process MC Particle
 706          fMCStack = fMCEvent->Stack();
 707          ProcessMCParticles();
 708       }
 709
 710
 711
 712       ProcessPhotonCandidates(); // Process this cuts gammas
 713       ProcessElectronCandidates(); // Process this cuts gammas
 714       CalculatePi0DalitzCandidates();
 715       CalculateBackground();
 716       UpdateEventByEventData();
 717
 718
 719
 720       fGoodGammas->Clear(); // delete this cuts good gammas
 721       fGoodVirtualGammas->Clear(); // delete this cuts good gammas
 722    }
 723
 724    fSelectorElectronIndex.clear();
 725    fSelectorPositronIndex.clear();
 726
 727    PostData( 1, fOutputContainer );
 728 }
 729
 730 void AliAnalysisTaskGammaConvDalitzV1::Terminate(const Option_t *)
 731 {
 732
 733    if( fElecSelector ){
 734
 735       if ( ((AliDalitzElectronCuts*)fElecSelector->GetDalitzElectronCuts())->GetCutHistograms() ){
 736          fOutputContainer->Add( ((AliDalitzElectronCuts*)fElecSelector->GetDalitzElectronCuts())->GetCutHistograms() );
 737       }
 738    }
 739
 740
 741    if ( fV0Reader ) {
 742
 743       if ( ((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms() ){
 744          fOutputContainer->Add(  ((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms() );
 745       }
 746    }
 747
 748    for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 749
 750       if( fCutElectronArray ){
 751          if( ((AliDalitzElectronCuts*)fCutElectronArray->At(iCut))->GetCutHistograms() ) {
 752             fCutFolder[iCut]->Add( ((AliDalitzElectronCuts*)fCutElectronArray->At(iCut))->GetCutHistograms() );
 753          }
 754       }
 755
 756       if( fCutMesonArray  ) {
 757          if( ((AliConversionMesonCuts*)fCutMesonArray->At(iCut))->GetCutHistograms() ) {
 758             fCutFolder[iCut]->Add( ((AliConversionMesonCuts*)fCutMesonArray->At(iCut))->GetCutHistograms());
 759          }
 760       }
 761
 762       if( fCutGammaArray ) {
 763          if( ((AliConversionCuts*)fCutGammaArray->At(iCut))->GetCutHistograms() ) {
 764             fCutFolder[iCut]->Add( ((AliConversionCuts*)fCutGammaArray->At(iCut))->GetCutHistograms()  );
 765          }
 766       }
 767    }
 768
 769 }
 770 //________________________________________________________________________
 771 void AliAnalysisTaskGammaConvDalitzV1::ProcessPhotonCandidates()
 772 {
 773    Int_t nV0 = 0;
 774    TList *GoodGammasStepOne = new TList();
 775    TList *GoodGammasStepTwo = new TList();
 776    // Loop over Photon Candidates allocated by ReaderV1
 777    for(Int_t i = 0; i < fReaderGammas->GetEntriesFast(); i++){
 778       AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) fReaderGammas->At(i);
 779       if(!PhotonCandidate) continue;
 780       if(!((AliConversionCuts*)fCutGammaArray->At(fiCut))->PhotonIsSelected(PhotonCandidate,fESDEvent)) continue;
 781
 782       if(!((AliConversionCuts*)fCutGammaArray->At(fiCut))->UseElecSharingCut() &&
 783          !((AliConversionCuts*)fCutGammaArray->At(fiCut))->UseToCloseV0sCut()){ // if no post reader loop is required add to events good gammas
 784          fGoodGammas->Add(PhotonCandidate);
 785          hESDConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
 786          if(fMCEvent){
 787             ProcessTruePhotonCandidates(PhotonCandidate);
 788          }
 789       }
 790       else if(((AliConversionCuts*)fCutGammaArray->At(fiCut))->UseElecSharingCut()){ // if Shared Electron cut is enabled, Fill array, add to step one
 791          ((AliConversionCuts*)fCutGammaArray->At(fiCut))->FillElectonLabelArray(PhotonCandidate,nV0);
 792          nV0++;
 793          GoodGammasStepOne->Add(PhotonCandidate);
 794       }
 795       else if(!((AliConversionCuts*)fCutGammaArray->At(fiCut))->UseElecSharingCut() &&
 796               ((AliConversionCuts*)fCutGammaArray->At(fiCut))->UseToCloseV0sCut()){ // shared electron is disabled, step one not needed -> step two
 797          GoodGammasStepTwo->Add(PhotonCandidate);
 798       }
 799    }
 800    if(((AliConversionCuts*)fCutGammaArray->At(fiCut))->UseElecSharingCut()){
 801       for(Int_t i = 0;i<GoodGammasStepOne->GetEntries();i++){
 802          AliAODConversionPhoton *PhotonCandidate= (AliAODConversionPhoton*) GoodGammasStepOne->At(i);
 803          if(!PhotonCandidate) continue;
 804          if(!((AliConversionCuts*)fCutGammaArray->At(fiCut))->RejectSharedElectronV0s(PhotonCandidate,i,GoodGammasStepOne->GetEntries())) continue;
 805          if(!((AliConversionCuts*)fCutGammaArray->At(fiCut))->UseToCloseV0sCut()){ // To Colse v0s cut diabled, step two not needed
 806             fGoodGammas->Add(PhotonCandidate);
 807             hESDConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
 808             if(fMCEvent){
 809                ProcessTruePhotonCandidates(PhotonCandidate);
 810             }
 811          }
 812          else GoodGammasStepTwo->Add(PhotonCandidate); // Close v0s cut enabled -> add to list two
 813       }
 814    }
 815    if(((AliConversionCuts*)fCutGammaArray->At(fiCut))->UseToCloseV0sCut()){
 816       for(Int_t i = 0;i<GoodGammasStepTwo->GetEntries();i++){
 817          AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) GoodGammasStepTwo->At(i);
 818          if(!PhotonCandidate) continue;
 819          if(!((AliConversionCuts*)fCutGammaArray->At(fiCut))->RejectToCloseV0s(PhotonCandidate,GoodGammasStepTwo,i)) continue;
 820          fGoodGammas->Add(PhotonCandidate); // Add gamma to current cut TList
 821          hESDConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt()); // Differences to old V0Reader in p_t due to conversion KF->TLorentzVector
 822          if(fMCEvent){
 823             ProcessTruePhotonCandidates(PhotonCandidate);
 824          }
 825       }
 826    }
 827
 828    delete GoodGammasStepOne;
 829    GoodGammasStepOne = 0x0;
 830    delete GoodGammasStepTwo;
 831    GoodGammasStepTwo = 0x0;
 832 }
 833
 834 //________________________________________________________________________
 835 void AliAnalysisTaskGammaConvDalitzV1::ProcessTruePhotonCandidates(AliAODConversionPhoton *TruePhotonCandidate)
 836 {
 837    // Process True Photons
 838    AliStack *MCStack = fMCEvent->Stack();
 839    TParticle *posDaughter = TruePhotonCandidate->GetPositiveMCDaughter(MCStack);
 840    TParticle *negDaughter = TruePhotonCandidate->GetNegativeMCDaughter(MCStack);
 841
 842    if(posDaughter == NULL || negDaughter == NULL) return; // One particle does not exist
 843    if(posDaughter->GetMother(0) != negDaughter->GetMother(0)){  // Not Same Mother == Combinatorial Bck
 844       return;
 845    }
 846    if(TMath::Abs(posDaughter->GetPdgCode())!=11 || TMath::Abs(negDaughter->GetPdgCode())!=11) return; //One Particle is not electron
 847    if(posDaughter->GetPdgCode()==negDaughter->GetPdgCode()) return; // Same Charge
 848    if(posDaughter->GetUniqueID() != 5 || negDaughter->GetUniqueID() !=5) return;// check if the daughters come from a conversion
 849
 850    TParticle *Photon = TruePhotonCandidate->GetMCParticle(MCStack);
 851    if(Photon->GetPdgCode() != 22) return; // Mother is no Photon
 852
 853    // True Photon
 854    hESDTrueConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
 855
 856    Int_t labelGamma = TruePhotonCandidate->GetMCParticleLabel(MCStack);
 857
 858     if( IsPi0DalitzDaughter(labelGamma) == kTRUE ) {
 859         hESDTruePi0DalitzConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
 860      }
 861
 862
 863
 864
 865 }
 866
 867
 868 //________________________________________________________________________
 869 void AliAnalysisTaskGammaConvDalitzV1::ProcessElectronCandidates(){
 870
 871    Double_t magField = fInputEvent->GetMagneticField();
 872
 873
 874    if( magField  < 0.0 ){
 875       magField =  1.0;
 876    }
 877    else {
 878       magField =  -1.0;
 879    }
 880
 881
 882
 883
 884    vector<Int_t> lGoodElectronIndex(0);
 885    vector<Int_t> lGoodPositronIndex(0);
 886
 887
 888
 889    for(UInt_t i = 0; i < fSelectorElectronIndex.size(); i++){
 890
 891       AliESDtrack* electronCandidate = fESDEvent->GetTrack(fSelectorElectronIndex[i]);
 892       if(! ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->ElectronIsSelected(electronCandidate) ) continue;
 893          lGoodElectronIndex.push_back(   fSelectorElectronIndex[i] );
 894          hESDDalitzElectronPt[fiCut]->Fill(electronCandidate->Pt());
 895             if( fMCEvent ) {
 896                     Int_t labelelectron = TMath::Abs( electronCandidate->GetLabel() );
 897                     if( labelelectron < fMCStack->GetNtrack() ){
 898                     TParticle* electron = fMCStack->Particle(labelelectron);
 899                       if( electron->GetPdgCode() ==  11 ){
 900                             hESDTrueElectronPt[fiCut]->Fill(electron->Pt());
 901                         if( IsPi0DalitzDaughter(labelelectron) == kTRUE ) {
 902                             hESDTruePi0DalitzElectronPt[fiCut]->Fill(electron->Pt());
 903                        }
 904                      }
 905                    }
 906             }
 907    }
 908
 909    for(UInt_t i = 0; i < fSelectorPositronIndex.size(); i++){
 910
 911       AliESDtrack* positronCandidate = fESDEvent->GetTrack( fSelectorPositronIndex[i] );
 912       if(! ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->ElectronIsSelected(positronCandidate) ) continue;
 913         lGoodPositronIndex.push_back(   fSelectorPositronIndex[i] );
 914         hESDDalitzPositronPt[fiCut]->Fill(positronCandidate->Pt());
 915         if( fMCEvent ) {
 916                     Int_t labelpositron = TMath::Abs( positronCandidate->GetLabel() );
 917                     if( labelpositron < fMCStack->GetNtrack() ) {
 918                     TParticle* positron = fMCStack->Particle(labelpositron);
 919                       if( positron->GetPdgCode() ==  -11 ){
 920                             hESDTruePositronPt[fiCut]->Fill(positron->Pt());
 921                         if( IsPi0DalitzDaughter(labelpositron) == kTRUE ) {
 922                             hESDTruePi0DalitzPositronPt[fiCut]->Fill(positron->Pt());
 923                        }
 924                      }
 925                }
 926         }
 927
 928    }
 929
 930
 931    vector<Bool_t> lElectronPsiIndex(lGoodElectronIndex.size(), kTRUE);
 932    vector<Bool_t> lPositronPsiIndex(lGoodPositronIndex.size(), kTRUE);
 933
 934
 935    if( ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->DoPsiPairCut() == kTRUE ){
 936
 937       for( UInt_t i = 0; i < lGoodElectronIndex.size(); i++ ) {
 938
 939          AliESDtrack *electronCandidate = fESDEvent->GetTrack(lGoodElectronIndex[i]);
 940
 941          for(UInt_t j = 0; j <  lGoodPositronIndex.size(); j++){
 942             AliESDtrack *positronCandidate = fESDEvent->GetTrack(lGoodPositronIndex[j]);
 943             Double_t psiPair = GetPsiPair(positronCandidate,electronCandidate);
 944             Double_t deltaPhi = magField * TVector2::Phi_mpi_pi( electronCandidate->GetConstrainedParam()->Phi()-positronCandidate->GetConstrainedParam()->Phi());
 945
 946             if( ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->IsFromGammaConversion(psiPair,deltaPhi ) ){
 947                lElectronPsiIndex[i] = kFALSE;
 948                lPositronPsiIndex[j] = kFALSE;
 949             }
 950          }
 951       }
 952    }
 953
 954
 955
 956
 957
 958
 959    for(UInt_t i = 0; i < lGoodElectronIndex.size(); i++){
 960
 961       if( lElectronPsiIndex[i] == kFALSE ) continue;
 962
 963
 964       AliESDtrack *electronCandidate = fESDEvent->GetTrack(lGoodElectronIndex[i]);
 965
 966       AliKFParticle electronCandidateKF( *electronCandidate->GetConstrainedParam(), ::kElectron );
 967
 968       for(UInt_t j = 0; j < lGoodPositronIndex.size(); j++){
 969
 970          if( lPositronPsiIndex[j] == kFALSE ) continue;
 971
 972          AliESDtrack *positronCandidate = fESDEvent->GetTrack(lGoodPositronIndex[j]);
 973          AliKFParticle positronCandidateKF( *positronCandidate->GetConstrainedParam(), ::kPositron );
 974          Bool_t isPhoton    = kFALSE;
 975          Bool_t isPi0Dalitz = kFALSE;
 976          Bool_t isEtaDalitz = kFALSE;
 977
 978
 979          Double_t psiPair = GetPsiPair(positronCandidate,electronCandidate);
 980          Double_t deltaPhi = magField * TVector2::Phi_mpi_pi( electronCandidate->GetConstrainedParam()->Phi()-positronCandidate->GetConstrainedParam()->Phi());
 981
 982
 983
 984          AliKFConversionPhoton* virtualPhoton = new AliKFConversionPhoton(electronCandidateKF,positronCandidateKF);
 985
 986
 987          AliKFVertex primaryVertexImproved(*fInputEvent->GetPrimaryVertex());
 988          primaryVertexImproved+=*virtualPhoton;
 989          virtualPhoton->SetProductionVertex(primaryVertexImproved);
 990
 991          virtualPhoton->SetTrackLabels( lGoodPositronIndex[j], lGoodElectronIndex[i]);
 992
 993          if(fMCEvent){
 994
 995             // AliStack *fMCStack= fMCEvent->Stack();
 996             Int_t labeln=TMath::Abs(electronCandidate->GetLabel());
 997             Int_t labelp=TMath::Abs(positronCandidate->GetLabel());
 998             TParticle *fNegativeMCParticle = fMCStack->Particle(labeln);
 999             TParticle *fPositiveMCParticle = fMCStack->Particle(labelp);
1000             if( fPositiveMCParticle && fNegativeMCParticle) {
1001                virtualPhoton->SetMCLabelPositive(labelp);
1002                virtualPhoton->SetMCLabelNegative(labeln);
1003
1004             }
1005
1006             TParticle *mcVgamma=virtualPhoton->GetMCParticle(fMCStack);
1007             Int_t labelgamma    = -1;
1008             Int_t labelelectron = -1;
1009             Int_t labelpositron = -1;
1010
1011             if(mcVgamma){
1012               // Check if it is a true photon
1013               if(mcVgamma->GetPdgCode() == 22){
1014                 isPhoton = kTRUE;
1015               }
1016               else if( IsDalitz( mcVgamma,labelgamma,labelelectron,labelpositron ) ){
1017                   if     ( mcVgamma->GetPdgCode() == 111 ) isPi0Dalitz = kTRUE;
1018                   else if( mcVgamma->GetPdgCode() == 221 ) isEtaDalitz = kTRUE;
1019               }
1020             }
1021            }
1022
1023
1024          AliAODConversionPhoton *vphoton = new AliAODConversionPhoton(virtualPhoton); //To Apply PsiPairCut
1025
1026          hESDEposEnegPsiPairDPhi[fiCut]->Fill(deltaPhi,psiPair);
1027          hESDEposEnegInvMassPt[fiCut]->Fill(vphoton->GetMass(),vphoton->Pt());
1028
1029          if( fMCEvent ) {
1030             if(isPhoton)        hESDEposEnegTruePhotonInvMassPt[fiCut]->Fill(vphoton->GetMass(),vphoton->Pt());
1031             else if(isPi0Dalitz)hESDEposEnegTruePi0DalitzInvMassPt[fiCut]->Fill(vphoton->GetMass(),vphoton->Pt());
1032             else if(isEtaDalitz)hESDEposEnegTrueEtaDalitzInvMassPt[fiCut]->Fill(vphoton->GetMass(),vphoton->Pt());
1033
1034          }
1035
1036          fGoodVirtualGammas->Add(  vphoton );
1037
1038       }
1039    }
1040
1041
1042    //Computing mixing event
1043
1044    for(UInt_t i = 0; i < lGoodElectronIndex.size(); i++){
1045
1046          if( lElectronPsiIndex[i] == kFALSE ) continue;
1047
1048           AliESDtrack *electronCandidate1 = fESDEvent->GetTrack(lGoodElectronIndex[i]);
1049
1050           AliKFParticle electronCandidate1KF( *electronCandidate1->GetConstrainedParam(), ::kElectron );
1051
1052
1053           for(UInt_t j = i+1; j < lGoodElectronIndex.size(); j++){
1054
1055                if( lElectronPsiIndex[j] == kFALSE ) continue;
1056
1057
1058                 AliESDtrack *electronCandidate2 = fESDEvent->GetTrack(lGoodElectronIndex[j]);
1059
1060                 AliKFParticle electronCandidate2KF( *electronCandidate2->GetConstrainedParam(), ::kElectron );
1061
1062                 AliKFConversionPhoton* virtualPhoton = new AliKFConversionPhoton(electronCandidate1KF,electronCandidate2KF);
1063
1064                 AliKFVertex primaryVertexImproved(*fInputEvent->GetPrimaryVertex());
1065                 primaryVertexImproved+=*virtualPhoton;
1066                 virtualPhoton->SetProductionVertex(primaryVertexImproved);
1067
1068
1069                 AliAODConversionPhoton *vphoton = new AliAODConversionPhoton(virtualPhoton);
1070                 hESDEposEnegLikeSignBackInvMassPt[fiCut]->Fill(vphoton->GetMass(),vphoton->Pt());
1071                 delete vphoton;
1072                 delete virtualPhoton;
1073
1074
1075           }
1076    }
1077
1078
1079    for(UInt_t i = 0; i < lGoodPositronIndex.size(); i++){
1080
1081          if( lPositronPsiIndex[i] == kFALSE ) continue;
1082
1083           AliESDtrack *positronCandidate1 = fESDEvent->GetTrack(lGoodPositronIndex[i]);
1084
1085           AliKFParticle positronCandidate1KF( *positronCandidate1->GetConstrainedParam(), ::kPositron );
1086
1087
1088           for(UInt_t j = i+1; j < lGoodPositronIndex.size(); j++){
1089
1090                if( lPositronPsiIndex[j] == kFALSE ) continue;
1091
1092                 AliESDtrack *positronCandidate2 = fESDEvent->GetTrack(lGoodPositronIndex[j]);
1093
1094                 AliKFParticle positronCandidate2KF( *positronCandidate2->GetConstrainedParam(), ::kPositron );
1095
1096                 AliKFConversionPhoton* virtualPhoton = new AliKFConversionPhoton(positronCandidate1KF,positronCandidate2KF);
1097                 AliKFVertex primaryVertexImproved(*fInputEvent->GetPrimaryVertex());
1098                 primaryVertexImproved+=*virtualPhoton;
1099                 virtualPhoton->SetProductionVertex(primaryVertexImproved);
1100
1101                 AliAODConversionPhoton *vphoton = new AliAODConversionPhoton(virtualPhoton);
1102                 hESDEposEnegLikeSignBackInvMassPt[fiCut]->Fill(vphoton->GetMass(),vphoton->Pt());
1103
1104
1105                 delete vphoton;
1106                 delete virtualPhoton;
1107
1108           }
1109    }
1110
1111
1112 }
1113
1114 //________________________________________________________________________
1115 void AliAnalysisTaskGammaConvDalitzV1::CalculatePi0DalitzCandidates(){
1116
1117    // Conversion Gammas
1118
1119
1120
1121
1122    if( fGoodGammas->GetEntries() > 0 && fGoodVirtualGammas->GetEntries() > 0 ){
1123
1124
1125       for(Int_t GammaIndex=0; GammaIndex<fGoodGammas->GetEntries(); GammaIndex++){
1126
1127          AliAODConversionPhoton *gamma=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(GammaIndex));
1128          if (gamma==NULL) continue;
1129          for(Int_t virtualGammaIndex=0;virtualGammaIndex<fGoodVirtualGammas->GetEntries();virtualGammaIndex++){
1130
1131             AliAODConversionPhoton *Vgamma=dynamic_cast<AliAODConversionPhoton*>(fGoodVirtualGammas->At(virtualGammaIndex));
1132             if (Vgamma==NULL) continue;
1133             //Check for same Electron ID
1134             if(gamma->GetTrackLabelPositive() == Vgamma->GetTrackLabelPositive() ||
1135                gamma->GetTrackLabelNegative() == Vgamma->GetTrackLabelNegative() ||
1136                gamma->GetTrackLabelNegative() == Vgamma->GetTrackLabelPositive() ||
1137                gamma->GetTrackLabelPositive() == Vgamma->GetTrackLabelNegative() ) continue;
1138
1139             AliAODConversionMother *pi0cand = new AliAODConversionMother(gamma,Vgamma);
1140             pi0cand->SetLabels(GammaIndex,virtualGammaIndex);
1141
1142             if( ( ((AliConversionMesonCuts*)fCutMesonArray->At(fiCut))->MesonIsSelected(pi0cand,kTRUE)) ){
1143                hESDMotherInvMassPt[fiCut]->Fill(pi0cand->M(),pi0cand->Pt());
1144                hESDPi0MotherInvMassPt[fiCut]->Fill(pi0cand->M(),pi0cand->Pt());
1145
1146                Double_t diffMass = pi0cand->M() - Vgamma->GetMass();
1147
1148                hESDPi0MotherDiffInvMassPt[fiCut]->Fill( diffMass , pi0cand->Pt() );
1149                // if(pi0cand->GetAlpha()<0.1){
1150                //    hESDMotherInvMassEalpha[fiCut]->Fill(pi0cand->M(),pi0cand->E());
1151                // }
1152                Int_t zbin= fBGHandler[fiCut]->GetZBinIndex(fESDEvent->GetPrimaryVertex()->GetZ());
1153                Int_t mbin = 0;
1154                if( ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->UseTrackMultiplicity() ){
1155                   mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fNumberOfESDTracks);
1156                } else {
1157                   mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGoodGammas->GetEntries());
1158                }
1159                Double_t sparesFill[4] = {pi0cand->M(),pi0cand->Pt(),(Double_t)zbin,(Double_t)mbin};
1160                sESDMotherInvMassPtZM[fiCut]->Fill(sparesFill,1);
1161                if(fMCEvent){
1162                   ProcessTrueMesonCandidates(pi0cand,gamma,Vgamma);
1163                }
1164             }
1165             delete pi0cand;
1166             pi0cand=0x0;
1167          }
1168       }
1169    }
1170 }
1171
1172 //________________________________________________________________________
1173 void AliAnalysisTaskGammaConvDalitzV1::CalculateBackground(){
1174
1175    Int_t zbin= fBGHandler[fiCut]->GetZBinIndex(fESDEvent->GetPrimaryVertex()->GetZ());
1176    Int_t mbin = 0;
1177
1178    Int_t method = 0;
1179
1180    method = ((AliDalitzElectronCuts*) fCutElectronArray->At(fiCut))->GetBKGMethod();
1181
1182
1183    if(((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->UseTrackMultiplicity()){
1184       mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fNumberOfESDTracks);
1185    } else {
1186       mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGoodGammas->GetEntries());
1187    }
1188
1189    if( method == 1 || method == 2 ) {
1190
1191       AliGammaConversionAODBGHandler::GammaConversionVertex *bgEventVertex = NULL;
1192
1193       if( ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->UseTrackMultiplicity() ) {
1194
1195          for(Int_t nEventsInBG=0;nEventsInBG<fBGHandler[fiCut]->GetNBGEvents();nEventsInBG++){
1196
1197             AliGammaConversionAODVector *previousEventV0s = fBGHandler[fiCut]->GetBGGoodV0s(zbin,mbin,nEventsInBG);
1198
1199             if(fMoveParticleAccordingToVertex == kTRUE && method == 1){
1200                bgEventVertex = fBGHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG);
1201             }
1202
1203             for(Int_t iCurrent=0;iCurrent<fGoodVirtualGammas->GetEntries();iCurrent++){
1204                AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGoodVirtualGammas->At(iCurrent));
1205
1206                for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
1207                   AliAODConversionPhoton previousGoodV0 = (AliAODConversionPhoton)(*(previousEventV0s->at(iPrevious)));
1208
1209                   if(fMoveParticleAccordingToVertex == kTRUE && method == 1 ){
1210                      MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
1211                   }
1212
1213                   AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&previousGoodV0);
1214                   if( ( ((AliConversionMesonCuts*)fCutMesonArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE) ) ){
1215                      hESDMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
1216                      Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin};
1217                      sESDMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
1218                   }
1219                   delete backgroundCandidate;
1220                   backgroundCandidate = 0x0;
1221                }
1222             }
1223          }
1224       }
1225       else{
1226          for(Int_t nEventsInBG=0;nEventsInBG <fBGHandler[fiCut]->GetNBGEvents();nEventsInBG++){
1227             AliGammaConversionAODVector *previousEventV0s = fBGHandler[fiCut]->GetBGGoodV0s(zbin,mbin,nEventsInBG);
1228             if(previousEventV0s){
1229                if(fMoveParticleAccordingToVertex == kTRUE && method == 1){
1230                   bgEventVertex = fBGHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG);
1231                }
1232                for(Int_t iCurrent=0;iCurrent<fGoodVirtualGammas->GetEntries();iCurrent++){
1233
1234                   AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGoodVirtualGammas->At(iCurrent));
1235
1236                   for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
1237
1238                      AliAODConversionPhoton previousGoodV0 = (AliAODConversionPhoton)(*(previousEventV0s->at(iPrevious)));
1239
1240                      if(fMoveParticleAccordingToVertex == kTRUE && method ==1){
1241                         MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
1242                      }
1243
1244                      AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&previousGoodV0);
1245                      if((((AliConversionMesonCuts*)fCutMesonArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE))){
1246                         hESDMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
1247                         Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin};
1248                         sESDMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
1249                      }
1250                      delete backgroundCandidate;
1251                      backgroundCandidate = 0x0;
1252                   }
1253                }
1254             }
1255          }
1256       }
1257    }
1258
1259    else if( method == 3 ){
1260
1261       for(Int_t iCurrent=0;iCurrent<fGoodVirtualGammas->GetEntries();iCurrent++){
1262
1263          AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGoodVirtualGammas->At(iCurrent));
1264
1265          for(Int_t iPrevious=0;iPrevious<fGammasPool[fiCut]->GetEntries();iPrevious++){
1266
1267             AliAODConversionPhoton previousGoodV0 = *(AliAODConversionPhoton*)((fGammasPool[fiCut]->At(iPrevious) ));
1268
1269
1270             AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&previousGoodV0);
1271
1272             if((((AliConversionMesonCuts*)fCutMesonArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE))){
1273
1274                hESDMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
1275                Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin};
1276                sESDMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
1277             }
1278             delete backgroundCandidate;
1279             backgroundCandidate = 0x0;
1280          }
1281       }
1282    }
1283
1284 }
1285 //________________________________________________________________________
1286 void AliAnalysisTaskGammaConvDalitzV1::UpdateEventByEventData(){
1287    //see header file for documentation
1288
1289    Int_t method = 0;
1290
1291    method = ( (AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->GetBKGMethod();
1292
1293    if( method == 1 ) {
1294
1295       if(fGoodGammas->GetEntries() > 0 ){
1296
1297          if( ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->UseTrackMultiplicity() ){
1298             fBGHandler[fiCut]->AddEvent(fGoodGammas,fESDEvent->GetPrimaryVertex()->GetX(),fESDEvent->GetPrimaryVertex()->GetY(),fESDEvent->GetPrimaryVertex()->GetZ(),fNumberOfESDTracks);
1299          }
1300
1301          else{ // means we use #V0s for multiplicity
1302             fBGHandler[fiCut]->AddEvent(fGoodGammas,fESDEvent->GetPrimaryVertex()->GetX(),fESDEvent->GetPrimaryVertex()->GetY(),fESDEvent->GetPrimaryVertex()->GetZ(),fGoodGammas->GetEntries());
1303          }
1304       }
1305    }
1306
1307    else if ( method == 2 ){
1308
1309       if(fGoodVirtualGammas->GetEntries() > 0 ){
1310          if(((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->UseTrackMultiplicity()){
1311             fBGHandler[fiCut]->AddEvent(fGoodVirtualGammas,fESDEvent->GetPrimaryVertex()->GetX(),fESDEvent->GetPrimaryVertex()->GetY(),fESDEvent->GetPrimaryVertex()->GetZ(),fNumberOfESDTracks);
1312          }
1313          else{ // means we use #V0s for multiplicity
1314             fBGHandler[fiCut]->AddEvent(fGoodVirtualGammas,fESDEvent->GetPrimaryVertex()->GetX(),fESDEvent->GetPrimaryVertex()->GetY(),fESDEvent->GetPrimaryVertex()->GetZ(),fGoodVirtualGammas->GetEntries());
1315          }
1316       }
1317    }
1318    else if ( method  == 3 ) {
1319
1320
1321
1322       for(Int_t index = 0; index < fGoodGammas->GetEntries(); index++){
1323
1324
1325          if ( fGammasPool[fiCut]->GetEntries() > ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->NumberOfRotationEvents() ){
1326             fGammasPool[fiCut]->RemoveLast();
1327          }
1328          fGammasPool[fiCut]->AddFirst( new AliAODConversionPhoton(*(AliAODConversionPhoton*)(fGoodGammas->At(index)) ) );
1329
1330       }
1331    }
1332 }
1333 //______________________________________________________________________
1334 void AliAnalysisTaskGammaConvDalitzV1::ProcessTrueMesonCandidates(AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate, AliAODConversionPhoton *TrueVirtualGammaCandidate)
1335 {
1336
1337    // Process True Mesons
1338
1339    AliStack *MCStack = fMCEvent->Stack();
1340
1341
1342    if(TrueGammaCandidate->GetV0Index()<fESDEvent->GetNumberOfV0s()){
1343
1344
1345       Bool_t isTruePi0 = kFALSE;
1346       Bool_t isTrueEta = kFALSE;
1347       Int_t gammaMCLabel = TrueGammaCandidate->GetMCParticleLabel(MCStack);
1348       Int_t gammaMotherLabel = -1;
1349
1350
1351
1352
1353       if(gammaMCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1354
1355
1356          // Daughters Gamma 0
1357          TParticle * negativeMC = (TParticle*)TrueGammaCandidate->GetNegativeMCDaughter(MCStack);
1358          TParticle * positiveMC = (TParticle*)TrueGammaCandidate->GetPositiveMCDaughter(MCStack);
1359          TParticle * gammaMC = (TParticle*)MCStack->Particle(gammaMCLabel);
1360
1361
1362          if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1363
1364             if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
1365
1366                if(gammaMC->GetPdgCode() == 22){ // ... with Gamma Mother
1367                   gammaMotherLabel=gammaMC->GetFirstMother();
1368                }
1369             }
1370          }
1371       }
1372
1373
1374       Int_t virtualGammaMCLabel = TrueVirtualGammaCandidate->GetMCParticleLabel(MCStack);
1375       Int_t virtualGammaMotherLabel = -1;
1376       Int_t virtualGamma = 1;
1377
1378       if(virtualGammaMCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1379          // Daughters Gamma 1
1380          TParticle * negativeMC = (TParticle*)TrueVirtualGammaCandidate->GetNegativeMCDaughter(MCStack);
1381          TParticle * positiveMC = (TParticle*)TrueVirtualGammaCandidate->GetPositiveMCDaughter(MCStack);
1382          TParticle * virtualGammaMotherMC = (TParticle*)MCStack->Particle(virtualGammaMCLabel);
1383
1384          if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1385
1386             if( virtualGammaMotherMC->GetPdgCode() != 22 ){
1387                virtualGammaMotherLabel=virtualGammaMCLabel;
1388             }
1389
1390             else if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
1391
1392                virtualGammaMotherLabel=virtualGammaMotherMC->GetFirstMother();
1393                virtualGamma = 0; //no virtual gamma
1394
1395             }
1396          }
1397       }
1398
1399
1400       if(gammaMotherLabel>=0 && ( gammaMotherLabel == virtualGammaMotherLabel) ){
1401
1402          if(((TParticle*)MCStack->Particle(virtualGammaMotherLabel))->GetPdgCode() == 111){
1403             isTruePi0=kTRUE;
1404          }
1405
1406          if(((TParticle*)MCStack->Particle(virtualGammaMotherLabel))->GetPdgCode() == 221){
1407             isTrueEta=kTRUE;
1408          }
1409
1410       }
1411
1412       if(isTruePi0 || isTrueEta ){ // True Pion or Eta
1413          if ( virtualGamma == 1 ) { //True Dalitz
1414             hESDTrueMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1415             if(gammaMotherLabel <= MCStack->GetNprimary()){ // Only primary pi0 for efficiency calculation
1416                hESDTruePrimaryMotherInvMassMCPt[fiCut]->Fill(Pi0Candidate->M(),((TParticle*)MCStack->Particle(virtualGammaMotherLabel))->Pt());
1417                if(isTruePi0){ // Only primaries for unfolding
1418                   hESDTruePrimaryPi0DalitzESDPtMCPt[fiCut]->Fill(Pi0Candidate->Pt(),((TParticle*)MCStack->Particle(virtualGammaMotherLabel))->Pt());
1419                }
1420             }
1421             if(gammaMotherLabel > MCStack->GetNprimary()){ // Secondary Meson
1422                hESDTrueSecondaryMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1423                //if (((TParticle*)MCStack->Particle(virtualGammaMotherLabel))->GetMother(0) >-1){
1424                //   if(MCStack->Particle(((TParticle*)MCStack->Particle(virtualGammaMotherLabel))->GetMother(0))->GetPdgCode()==kK0Short){
1425                //      hESDTrueSecondaryMotherFromK0sInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1426                //   }
1427             }
1428          }
1429          else if ( virtualGamma == 0 ){
1430             hESDTrueMotherPi0GGInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt()); // Pi0 from GG
1431          }
1432       }
1433
1434       if(!isTruePi0 && !isTrueEta){ // Background
1435          if(gammaMotherLabel>-1 && virtualGammaMotherLabel>-1 && virtualGamma == 0){ // Both Tracks are Photons and have a mother but not Pi0 or Eta
1436             hESDTrueBckGGInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1437          } else { // No photon or without mother
1438             hESDTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1439          }
1440       }
1441    }
1442 }
1443 //________________________________________________________________________
1444 void AliAnalysisTaskGammaConvDalitzV1::MoveParticleAccordingToVertex(AliAODConversionPhoton* particle,const AliGammaConversionAODBGHandler::GammaConversionVertex *vertex){
1445    //see header file for documentation
1446
1447    Double_t dx = vertex->fX - fESDEvent->GetPrimaryVertex()->GetX();
1448    Double_t dy = vertex->fY - fESDEvent->GetPrimaryVertex()->GetY();
1449    Double_t dz = vertex->fZ - fESDEvent->GetPrimaryVertex()->GetZ();
1450
1451    Double_t movedPlace[3] = {particle->GetConversionX() - dx,particle->GetConversionY() - dy,particle->GetConversionZ() - dz};
1452    particle->SetConversionPoint(movedPlace);
1453 }
1454
1455
1456 //________________________________________________________________________
1457 void AliAnalysisTaskGammaConvDalitzV1::CountESDTracks(){
1458
1459    // Using standard function for setting Cuts
1460    Bool_t selectPrimaries=kTRUE;
1461    AliESDtrackCuts *EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
1462    EsdTrackCuts->SetMaxDCAToVertexZ(2);
1463    EsdTrackCuts->SetEtaRange(-0.8, 0.8);
1464    EsdTrackCuts->SetPtRange(0.15);
1465
1466    fNumberOfESDTracks = 0;
1467    for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
1468       AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
1469       if(!curTrack) continue;
1470       if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
1471    }
1472    delete EsdTrackCuts;
1473    EsdTrackCuts=0x0;
1474
1475    return;
1476 }
1477
1478 //_____________________________________________________________________________
1479 void AliAnalysisTaskGammaConvDalitzV1::ProcessMCParticles()
1480 {
1481
1482    // Loop over all primary MC particle
1483
1484    for(Int_t i = 0; i < fMCStack->GetNprimary(); i++) {
1485       TParticle* particle = (TParticle *)fMCStack->Particle(i);
1486       if (!particle) continue;
1487
1488
1489     if(((AliConversionCuts*)fCutGammaArray->At(fiCut))->PhotonIsSelectedMC(particle,fMCStack,kFALSE)){
1490          hMCAllGammaPt[fiCut]->Fill(particle->Pt()); // All MC Gamma
1491     }
1492
1493     if(  ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->ElectronIsSelectedMC(i,fMCStack)){
1494         if( particle->GetPdgCode() == -11)hMCAllPositronsPt[fiCut]->Fill(particle->Pt()); // All positrons
1495         if( particle->GetPdgCode() ==  11)hMCAllElectronsPt[fiCut]->Fill(particle->Pt()); // All electrons
1496     }
1497
1498     if(((AliConversionMesonCuts*)fCutMesonArray->At(fiCut))->MesonIsSelectedMC(particle,fMCStack)){
1499
1500             if(particle->GetPdgCode() == 111)hMCPi0GGPt[fiCut]->Fill( particle->Pt() ); // All MC Pi0 GG decay
1501             if(particle->GetPdgCode() == 221)hMCEtaGGPt[fiCut]->Fill( particle->Pt() ); // All MC Eta GG decay
1502     }
1503
1504     if(((AliConversionMesonCuts*)fCutMesonArray->At(fiCut))->MesonIsSelectedMCDalitz(particle,fMCStack)){
1505
1506          if(particle->GetPdgCode() == 111)hMCPi0Pt[fiCut]->Fill(particle->Pt()); // All MC Pi0
1507          if(particle->GetPdgCode() == 221)hMCEtaPt[fiCut]->Fill(particle->Pt()); // All MC Eta
1508
1509          // Check the acceptance for gamma and electrons
1510
1511          Int_t labelgamma =    -1;
1512          Int_t labelelectron = -1;
1513          Int_t labelpositron = -1;
1514
1515         if( IsDalitz( particle,labelgamma,labelelectron,labelpositron) == kTRUE ) {
1516
1517
1518                TParticle *gamma    = fMCStack->Particle(labelgamma);
1519                TParticle *electron = fMCStack->Particle(labelelectron);
1520                TParticle *positron = fMCStack->Particle(labelpositron);
1521
1522
1523                if(((AliConversionCuts*)fCutGammaArray->At(fiCut))->PhotonIsSelectedMC(gamma,fMCStack,kFALSE) &&
1524                   TMath::Abs( electron->Eta() ) < ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->GetEtaCut()  &&
1525                   TMath::Abs( positron->Eta() ) < ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->GetEtaCut() ){
1526                   if(particle->GetPdgCode() == 111){
1527                         hMCPi0InAccPt[fiCut]->Fill(particle->Pt()); // MC Pi0Dalitz with gamma and e+e- in acc
1528                         hMCPi0DalitzGammaPt[fiCut]->Fill( gamma->Pt());
1529                         hMCPi0DalitzPositronPt[fiCut]->Fill( positron->Pt());
1530                         hMCPi0DalitzElectronPt[fiCut]->Fill( electron->Pt());
1531                   }
1532                   if(particle->GetPdgCode() == 221)hMCEtaInAccPt[fiCut]->Fill(particle->Pt()); // MC EtaDalitz with gamma and e+e- in acc
1533                }
1534          }
1535       }
1536    }
1537 }
1538 //_____________________________________________________________________________
1539 Bool_t AliAnalysisTaskGammaConvDalitzV1::IsDalitz(TParticle *fMCMother,Int_t &labelgamma, Int_t &labelelectron,Int_t &labelpositron)
1540 {
1541
1542             if( fMCMother->GetNDaughters() != 3 ) return kFALSE;
1543             if( fMCMother->GetPdgCode() != 111 && fMCMother->GetPdgCode() != 221 ) return kFALSE;
1544
1545             TParticle *positron = 0x0;
1546             TParticle *electron = 0x0;
1547             TParticle *gamma    = 0x0;
1548
1549
1550
1551
1552             for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
1553
1554
1555                TParticle* temp = (TParticle*)fMCStack->Particle( index );
1556
1557                switch( temp->GetPdgCode() ) {
1558                case ::kPositron:
1559                                 positron =  temp;
1560                                 labelpositron = index;
1561                                 break;
1562                case ::kElectron:
1563                                 electron =  temp;
1564                                 labelelectron = index;
1565                                 break;
1566                case ::kGamma:
1567                                 gamma    =  temp;
1568                                 labelgamma = index;
1569                                 break;
1570                }
1571             }
1572
1573             if( positron && electron && gamma) return kTRUE;
1574
1575 return kFALSE;
1576 }
1577 //_____________________________________________________________________________________
1578 Bool_t AliAnalysisTaskGammaConvDalitzV1::IsPi0DalitzDaughter( Int_t label ) const
1579 {
1580 //
1581 // Returns true if the particle comes from Pi0 -> e+ e- gamma
1582 //
1583         Bool_t ePlusFlag  = kFALSE;
1584         Bool_t eMinusFlag = kFALSE;
1585         Bool_t gammaFlag  = kFALSE;
1586
1587         Int_t motherLabel = fMCStack->Particle( label )->GetMother(0);
1588
1589         if( motherLabel < 0 ) return kFALSE;
1590
1591         TParticle* mother = fMCStack->Particle( motherLabel );
1592
1593         if ( mother->GetPdgCode() != ::kPi0 ) return kFALSE;
1594
1595         if ( mother->GetNDaughters() != 3 ) return kFALSE;
1596
1597         for( Int_t idx = mother->GetFirstDaughter(); idx <= mother->GetLastDaughter(); ++idx ) {
1598                 switch( fMCStack->Particle(idx)->GetPdgCode()) {
1599                         case ::kPositron:
1600                                 ePlusFlag  = kTRUE;
1601                                 break;
1602                         case ::kElectron:
1603                                 eMinusFlag = kTRUE;
1604                                 break;
1605                         case ::kGamma:
1606                                 gammaFlag  = kTRUE;
1607                                 break;
1608                 }
1609         }
1610         return ( ePlusFlag && eMinusFlag && gammaFlag );
1611 }
1612
1613
1614 //_____________________________________________________________________________
1615 Double_t AliAnalysisTaskGammaConvDalitzV1::GetPsiPair( const AliESDtrack *trackPos, const AliESDtrack *trackNeg ) const
1616 {
1617    //
1618    // This angle is a measure for the contribution of the opening in polar
1619    // direction Δ0 to the opening angle ξ Pair
1620    //
1621    // Ref. Measurement of photons via conversion pairs with the PHENIX experiment at RHIC
1622    //      Master Thesis. Thorsten Dahms. 2005
1623    // https://twiki.cern.ch/twiki/pub/ALICE/GammaPhysicsPublications/tdahms_thesis.pdf
1624    //
1625    Double_t momPos[3];
1626    Double_t momNeg[3];
1627    if( trackPos->GetConstrainedPxPyPz(momPos) == 0 ) trackPos->GetPxPyPz( momPos );
1628    if( trackNeg->GetConstrainedPxPyPz(momNeg) == 0 ) trackNeg->GetPxPyPz( momNeg );
1629
1630    TVector3 posDaughter;
1631    TVector3 negDaughter;
1632
1633    posDaughter.SetXYZ( momPos[0], momPos[1], momPos[2] );
1634    negDaughter.SetXYZ( momNeg[0], momNeg[1], momNeg[2] );
1635
1636    Double_t deltaTheta = negDaughter.Theta() - posDaughter.Theta();
1637    Double_t openingAngle =  posDaughter.Angle( negDaughter );  //TMath::ACos( posDaughter.Dot(negDaughter)/(negDaughter.Mag()*posDaughter.Mag()) );
1638
1639    if( openingAngle < 1e-20 ) return 0.;
1640
1641    Double_t psiAngle = TMath::ASin( deltaTheta/openingAngle );
1642
1643    return psiAngle;
1644 }