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
  29 #include "AliStack.h"
  30 #include "AliAnalysisManager.h"
  31 #include "AliESDInputHandler.h"
  32 #include "AliESDtrack.h"
  33 #include "AliMCEvent.h"
  34 #include "AliStack.h"
  35 #include "AliMCEventHandler.h"
  36 #include "AliPID.h"
  37 #include "AliLog.h"
  38 #include "AliESDtrackCuts.h"
  39 #include "AliESDpidCuts.h"
  40 #include "AliMCEvent.h"
  41 #include "AliESDv0.h"
  42 #include "AliESDEvent.h"
  43 #include "AliESDpid.h"
  44 #include "AliKFParticle.h"
  45 #include "AliMCEventHandler.h"
  46 #include "AliKFVertex.h"
  47 #include "AliTriggerAnalysis.h"
  48 #include "AliCentrality.h"
  49 #include "AliMultiplicity.h"
  50 #include "AliAnalysisTaskGammaConvDalitzV1.h"
  51 #include "TH1.h"
  52 #include "TH2F.h"
  53 #include "THnSparse.h"
  54
  55 ClassImp( AliAnalysisTaskGammaConvDalitzV1 )
  56
  57 //-----------------------------------------------------------------------------------------------
  58 AliAnalysisTaskGammaConvDalitzV1::AliAnalysisTaskGammaConvDalitzV1():
  59    fV0Reader(NULL),
  60    fElecSelector(NULL),
  61    fBGHandler(NULL),
  62    fESDEvent(NULL),
  63    fMCEvent(NULL),
  64    fMCStack(NULL),
  65    fCutFolder(NULL),
  66    fESDList(NULL),
  67    fBackList(NULL),
  68    fTrueList(NULL),
  69    fMCList(NULL),
  70    fOutputContainer(0),
  71    fReaderGammas(NULL),
  72    fSelectorElectronIndex(0),
  73    fSelectorPositronIndex(0),
  74    fGoodGammas(NULL),
  75    fGoodVirtualGammas(NULL),
  76    fGoodElectrons(NULL),
  77    fGoodPositrons(NULL),
  78    fCutGammaArray(NULL),
  79    fCutElectronArray(NULL),
  80    fCutMesonArray(NULL),
  81    fGammasPool(NULL),
  82    fConversionCuts(NULL),
  83    hESDConvGammaPt(NULL),
  84    hESDDalitzElectronPt(NULL),
  85    hESDDalitzPositronPt(NULL),
  86    hESDEposEnegPsiPairDPhi(NULL),
  87    hESDMotherInvMassPt(NULL),
  88    sESDMotherInvMassPtZM(NULL),
  89    hESDMotherBackInvMassPt(NULL),
  90    sESDMotherBackInvMassPtZM(NULL),
  91    hMCPi0Pt(NULL),
  92    hMCEtaPt(NULL),
  93    hMCPi0InAccPt(NULL),
  94    hMCEtaInAccPt(NULL),
  95    hESDTrueMotherInvMassPt(NULL),
  96    hESDTrueMotherPi0GGInvMassPt(NULL),
  97    hESDTruePrimaryMotherInvMassMCPt(NULL),
  98    hESDTruePrimaryPi0DalitzESDPtMCPt(NULL),
  99    hESDTrueSecondaryMotherInvMassPt(NULL),
 100    hESDTrueSecondaryMotherFromK0sInvMassPt(NULL),
 101    hESDTrueBckGGInvMassPt(NULL),
 102    hESDTrueBckContInvMassPt(NULL),
 103    hESDTrueMotherGGInvMassPt(NULL),
 104    hESDTrueConvGammaPt(NULL),
 105    hNEvents(NULL),
 106    hNGoodESDTracks(NULL),
 107    fRandom(0),
 108    fUnsmearedPx(NULL),
 109    fUnsmearedPy(NULL),
 110    fUnsmearedPz(NULL),
 111    fUnsmearedE(NULL),
 112    fnCuts(0),
 113    fiCut(0),
 114    fNumberOfESDTracks(0),
 115    fMoveParticleAccordingToVertex(kFALSE),
 116    fIsHeavyIon(kFALSE),
 117    fDoMesonAnalysis(kTRUE)
 118 {
 119    // Define input and output slots here
 120   // DefineInput(0, TChain::Class());
 121   // DefineOutput(1, TList::Class());
 122 }
 123
 124 //-----------------------------------------------------------------------------------------------
 125 AliAnalysisTaskGammaConvDalitzV1::AliAnalysisTaskGammaConvDalitzV1( const char* name ):
 126    AliAnalysisTaskSE(name),
 127    fV0Reader(NULL),
 128    fElecSelector(NULL),
 129    fBGHandler(NULL),
 130    fESDEvent(NULL),
 131    fMCEvent(NULL),
 132    fMCStack(NULL),
 133    fCutFolder(NULL),
 134    fESDList(NULL),
 135    fBackList(NULL),
 136    fTrueList(NULL),
 137    fMCList(NULL),
 138    fOutputContainer(0),
 139    fReaderGammas(NULL),
 140    fSelectorElectronIndex(0),
 141    fSelectorPositronIndex(0),
 142    fGoodGammas(NULL),
 143    fGoodVirtualGammas(NULL),
 144    fGoodElectrons(NULL),
 145    fGoodPositrons(NULL),
 146    fCutGammaArray(NULL),
 147    fCutElectronArray(NULL),
 148    fCutMesonArray(NULL),
 149    fGammasPool(NULL),
 150    fConversionCuts(NULL),
 151    hESDConvGammaPt(NULL),
 152    hESDDalitzElectronPt(NULL),
 153    hESDDalitzPositronPt(NULL),
 154    hESDEposEnegPsiPairDPhi(NULL),
 155    hESDMotherInvMassPt(NULL),
 156    sESDMotherInvMassPtZM(NULL),
 157    hESDMotherBackInvMassPt(NULL),
 158    sESDMotherBackInvMassPtZM(NULL),
 159    hMCPi0Pt(NULL),
 160    hMCEtaPt(NULL),
 161    hMCPi0InAccPt(NULL),
 162    hMCEtaInAccPt(NULL),
 163    hESDTrueMotherInvMassPt(NULL),
 164    hESDTrueMotherPi0GGInvMassPt(NULL),
 165    hESDTruePrimaryMotherInvMassMCPt(NULL),
 166    hESDTruePrimaryPi0DalitzESDPtMCPt(NULL),
 167    hESDTrueSecondaryMotherInvMassPt(NULL),
 168    hESDTrueSecondaryMotherFromK0sInvMassPt(NULL),
 169    hESDTrueBckGGInvMassPt(NULL),
 170    hESDTrueBckContInvMassPt(NULL),
 171    hESDTrueMotherGGInvMassPt(NULL),
 172    hESDTrueConvGammaPt(NULL),
 173    hNEvents(NULL),
 174    hNGoodESDTracks(NULL),
 175    fRandom(0),
 176    fUnsmearedPx(NULL),
 177    fUnsmearedPy(NULL),
 178    fUnsmearedPz(NULL),
 179    fUnsmearedE(NULL),
 180    fnCuts(0),
 181    fiCut(0),
 182    fNumberOfESDTracks(0),
 183    fMoveParticleAccordingToVertex(kFALSE),
 184    fIsHeavyIon(kFALSE),
 185    fDoMesonAnalysis(kTRUE)
 186 {
 187         DefineInput(0, TChain::Class());
 188         DefineOutput(1, TList::Class());
 189        // DefineOutput(2, TList::Class());
 190
 191 }
 192
 193 //-----------------------------------------------------------------------------------------------
 194 AliAnalysisTaskGammaConvDalitzV1::~AliAnalysisTaskGammaConvDalitzV1()
 195 {
 196 //
 197 // virtual destructor
 198 //
 199
 200    if(fGoodGammas){
 201       delete fGoodGammas;
 202       fGoodGammas = 0x0;
 203    }
 204    if(fGoodVirtualGammas){
 205       delete fGoodVirtualGammas;
 206       fGoodGammas = 0x0;
 207    }
 208    if(fGoodElectrons){
 209       delete fGoodGammas;
 210       fGoodGammas = 0x0;
 211    }
 212    if(fGoodPositrons){
 213       delete fGoodGammas;
 214       fGoodGammas = 0x0;
 215    }
 216    if(fBGHandler){
 217       delete[] fBGHandler;
 218       fBGHandler = 0x0;
 219    }
 220 }
 221
 222 //___________________________________________________________
 223 void AliAnalysisTaskGammaConvDalitzV1::InitBack(){
 224
 225    Double_t *zBinLimitsArray = new Double_t[9];
 226    zBinLimitsArray[0] = -50.00;
 227    zBinLimitsArray[1] = -3.375;
 228    zBinLimitsArray[2] = -1.605;
 229    zBinLimitsArray[3] = -0.225;
 230    zBinLimitsArray[4] = 1.065;
 231    zBinLimitsArray[5] = 2.445;
 232    zBinLimitsArray[6] = 4.245;
 233    zBinLimitsArray[7] = 50.00;
 234    zBinLimitsArray[8] = 1000.00;
 235
 236    Double_t *multiplicityBinLimitsArrayTracks = new Double_t[6];
 237    multiplicityBinLimitsArrayTracks[0] = 0;
 238    multiplicityBinLimitsArrayTracks[1] = 8.5;
 239    multiplicityBinLimitsArrayTracks[2] = 16.5;
 240    multiplicityBinLimitsArrayTracks[3] = 27.5;
 241    multiplicityBinLimitsArrayTracks[4] = 41.5;
 242    multiplicityBinLimitsArrayTracks[5] = 200.;
 243    if(fIsHeavyIon){
 244       multiplicityBinLimitsArrayTracks[0] = 0;
 245       multiplicityBinLimitsArrayTracks[1] = 200.;
 246       multiplicityBinLimitsArrayTracks[2] = 500.;
 247       multiplicityBinLimitsArrayTracks[3] = 1000.;
 248       multiplicityBinLimitsArrayTracks[4] = 1500.;
 249       multiplicityBinLimitsArrayTracks[5] = 5000.;
 250    }
 251
 252    Double_t *multiplicityBinLimitsArrayV0s = new Double_t[5];
 253    multiplicityBinLimitsArrayV0s[0] = 2;
 254    multiplicityBinLimitsArrayV0s[1] = 3;
 255    multiplicityBinLimitsArrayV0s[2] = 4;
 256    multiplicityBinLimitsArrayV0s[3] = 5;
 257    multiplicityBinLimitsArrayV0s[4] = 9999;
 258    if(fIsHeavyIon){
 259       multiplicityBinLimitsArrayV0s[0] = 2;
 260       multiplicityBinLimitsArrayV0s[1] = 10;
 261       multiplicityBinLimitsArrayV0s[2] = 30;
 262       multiplicityBinLimitsArrayV0s[3] = 50;
 263       multiplicityBinLimitsArrayV0s[4] = 9999;
 264    }
 265
 266    const Int_t nDim = 4;
 267    Int_t nBins[nDim] = {1000,250,8,5};
 268    Double_t xMin[nDim] = {0,0, 0,0};
 269    Double_t xMax[nDim] = {1,25,8,5};
 270
 271    sESDMotherInvMassPtZM = new THnSparseF*[fnCuts];
 272    sESDMotherBackInvMassPtZM = new THnSparseF*[fnCuts];
 273
 274    fBGHandler = new AliGammaConversionAODBGHandler*[fnCuts];
 275
 276    for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 277
 278
 279       TString cutstringElectron =((AliDalitzElectronCuts*)fCutElectronArray->At(iCut))->GetCutNumber();
 280       TString cutstringMeson= ((AliConversionMesonCuts*)fCutMesonArray->At(iCut))->GetCutNumber();
 281       TString cutstringGamma = ((AliConversionCuts*)fCutGammaArray->At(iCut))->GetCutNumber();
 282
 283
 284       fBackList[iCut] = new TList();
 285       fBackList[iCut]->SetName(Form("%s_%s_%s Back histograms",cutstringGamma.Data(),cutstringElectron.Data(),cutstringMeson.Data()));
 286       fBackList[iCut]->SetOwner(kTRUE);
 287       fCutFolder[iCut]->Add(fBackList[iCut]);
 288
 289
 290
 291
 292       sESDMotherInvMassPtZM[iCut] = new THnSparseF("Back_Mother_InvMass_Pt_z_m","Back_Mother_InvMass_Pt_z_m",nDim,nBins,xMin,xMax);
 293       sESDMotherInvMassPtZM[iCut]->Sumw2();
 294       fBackList[iCut]->Add(sESDMotherInvMassPtZM[iCut]);
 295       sESDMotherBackInvMassPtZM[iCut] = new THnSparseF("Back_Back_InvMass_Pt_z_m","Back_Back_InvMass_Pt_z_m",nDim,nBins,xMin,xMax);
 296       sESDMotherBackInvMassPtZM[iCut]->Sumw2();
 297       fBackList[iCut]->Add(sESDMotherBackInvMassPtZM[iCut]);
 298
 299       if(((AliConversionCuts*)fCutGammaArray->At(iCut))->UseTrackMultiplicity()){
 300          fBGHandler[iCut] = new AliGammaConversionAODBGHandler(9,6,((AliConversionCuts*)fCutGammaArray->At(iCut))->NumberOfRotationEvents());
 301          fBGHandler[iCut]->Initialize(zBinLimitsArray, multiplicityBinLimitsArrayTracks);
 302       }
 303       else{
 304          fBGHandler[iCut] = new AliGammaConversionAODBGHandler(9,5,((AliConversionCuts*)fCutGammaArray->At(iCut))->NumberOfRotationEvents());
 305          fBGHandler[iCut]->Initialize(zBinLimitsArray, multiplicityBinLimitsArrayV0s);
 306       }
 307    }
 308 }
 309
 310 //______________________________________________________________________
 311 void AliAnalysisTaskGammaConvDalitzV1::UserCreateOutputObjects()
 312 {
 313 //
 314 // Create ouput objects
 315 //
 316
 317         // Create the output container
 318         if(fOutputContainer != NULL){
 319                 delete fOutputContainer;
 320                 fOutputContainer = NULL;
 321         }
 322         if(fOutputContainer == NULL){
 323                 fOutputContainer = new TList();
 324                 fOutputContainer->SetOwner(kTRUE);
 325         }
 326
 327         fGoodGammas = new TList();
 328         //fGoodGammas->SetOwner(kTRUE);
 329
 330
 331         fGoodVirtualGammas = new TList();
 332            //fGoodVirtualGammas->SetOwner(kTRUE);
 333
 334
 335
 336
 337
 338         fGammasPool = new TList*[fnCuts];
 339
 340         fCutFolder = new TList*[fnCuts];
 341         fESDList = new TList*[fnCuts];
 342         fBackList = new TList*[fnCuts];
 343         hNEvents = new TH1I*[fnCuts];
 344         hNGoodESDTracks = new TH1I*[fnCuts];
 345         hESDConvGammaPt = new TH1F*[fnCuts];
 346         hESDDalitzElectronPt = new TH1F*[fnCuts];
 347         hESDDalitzPositronPt = new TH1F*[fnCuts];
 348         hESDEposEnegPsiPairDPhi = new TH2F*[fnCuts];
 349
 350         //if(fDoMesonAnalysis){
 351                 hESDMotherInvMassPt = new TH2F*[fnCuts];
 352                 hESDMotherBackInvMassPt = new TH2F*[fnCuts];
 353         //}
 354         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 355
 356
 357                 TString cutstringElectron =((AliDalitzElectronCuts*)fCutElectronArray->At(iCut))->GetCutNumber();
 358                 TString cutstringMeson= ((AliConversionMesonCuts*)fCutMesonArray->At(iCut))->GetCutNumber();
 359                 TString cutstringGamma = ((AliConversionCuts*)fCutGammaArray->At(iCut))->GetCutNumber();
 360
 361                 fCutFolder[iCut] = new TList();
 362                 fCutFolder[iCut]->SetName(Form("Cut Number %s_%s_%s",cutstringGamma.Data(),cutstringElectron.Data(),cutstringMeson.Data()));
 363                 fCutFolder[iCut]->SetOwner(kTRUE);
 364                 fOutputContainer->Add(fCutFolder[iCut]);
 365                 fESDList[iCut] = new TList();
 366                 fESDList[iCut]->SetName(Form("%s_%s_%s ESD histograms",cutstringGamma.Data(),cutstringElectron.Data(),cutstringMeson.Data()));
 367                 fESDList[iCut]->SetOwner(kTRUE);
 368                 hNEvents[iCut] = new TH1I("NEvents","NEvents",7,-0.5,6.5);
 369                 fESDList[iCut]->Add(hNEvents[iCut]);
 370                 if(fIsHeavyIon) hNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",3000,0,3000);
 371                 else hNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",200,0,200);
 372                 fESDList[iCut]->Add(hNGoodESDTracks[iCut]);
 373                 hESDConvGammaPt[iCut] = new TH1F("ESD_ConvGamma_Pt","ESD_ConvGamma_Pt",250,0,25);
 374                 fESDList[iCut]->Add(hESDConvGammaPt[iCut]);
 375                 hESDDalitzElectronPt[iCut] = new TH1F("ESD_DalitzElectron_Pt","ESD_DalitzElectron_Pt",250,0,25);
 376                 fESDList[iCut]->Add(hESDDalitzElectronPt[iCut]);
 377                 hESDDalitzPositronPt[iCut] = new TH1F("ESD_DalitzPositron_Pt","ESD_DalitzPositron_Pt",250,0,25);
 378                 fESDList[iCut]->Add(hESDDalitzPositronPt[iCut]);
 379
 380
 381                 hESDEposEnegPsiPairDPhi[iCut] = new TH2F("ESD_EposEneg_PsiPair_DPhi","ESD_EposEneg_PsiPair_DPhi", 100, -1.0*TMath::Pi(), TMath::Pi(),100,-1.0*TMath::Pi(), TMath::Pi() );
 382                 fESDList[iCut]->Add(hESDEposEnegPsiPairDPhi[iCut]);
 383
 384
 385
 386
 387                // hESDDalitzPositronPt[iCut] = new TH1F("ESD_EposEneg_PsiPair_vs_DPhi","ESD_EposEneg_PsiPair_vs_DPhi",250,0,25);
 388
 389 //              if(fDoMesonAnalysis){
 390                         hESDMotherInvMassPt[iCut] = new TH2F("ESD_DalitzMother_InvMass_Pt","ESD_DalitzMother_InvMass_Pt",1000,0,1,250,0,25);
 391                         fESDList[iCut]->Add(hESDMotherInvMassPt[iCut]);
 392                         hESDMotherBackInvMassPt[iCut] = new TH2F("ESD_DalitzBackground_InvMass_Pt","ESD_DalitzBackground_InvMass_Pt",1000,0,1,250,0,25);
 393                         fESDList[iCut]->Add(hESDMotherBackInvMassPt[iCut]);
 394 //              }
 395                 fCutFolder[iCut]->Add(fESDList[iCut]);
 396
 397
 398              if( fCutElectronArray ){
 399                 if( ((AliDalitzElectronCuts*)fCutElectronArray->At(iCut))->GetCutHistograms() ) {
 400                     fCutFolder[iCut]->Add( ((AliDalitzElectronCuts*)fCutElectronArray->At(iCut))->GetCutHistograms() );
 401                 }
 402             }
 403
 404              if( fCutMesonArray  ) {
 405                 if( ((AliConversionMesonCuts*)fCutMesonArray->At(iCut))->GetCutHistograms() ) {
 406                     fCutFolder[iCut]->Add( ((AliConversionMesonCuts*)fCutMesonArray->At(iCut))->GetCutHistograms());
 407                 }
 408              }
 409
 410             if( fCutGammaArray ) {
 411                 if( ((AliConversionCuts*)fCutGammaArray->At(iCut))->GetCutHistograms() ) {
 412                     fCutFolder[iCut]->Add( ((AliConversionCuts*)fCutGammaArray->At(iCut))->GetCutHistograms()  );
 413                 }
 414             }
 415
 416
 417         }
 418
 419         //if(fDoMesonAnalysis){
 420                 InitBack(); // Init Background Handler
 421         //}
 422
 423         if(MCEvent()){
 424                 // MC Histogramms
 425                 fMCList = new TList*[fnCuts];
 426                 // True Histogramms
 427                 fTrueList = new TList*[fnCuts];
 428                 hESDTrueConvGammaPt = new TH1F*[fnCuts];
 429                 //if(fDoMesonAnalysis){
 430                         hMCPi0Pt = new TH1F*[fnCuts];
 431                         hMCEtaPt = new TH1F*[fnCuts];
 432                         hMCPi0InAccPt = new TH1F*[fnCuts];
 433                         hMCEtaInAccPt = new TH1F*[fnCuts];
 434
 435                         hESDTrueMotherInvMassPt = new TH2F*[fnCuts];
 436                         hESDTrueMotherPi0GGInvMassPt = new TH2F*[fnCuts];
 437                         hESDTruePrimaryPi0DalitzESDPtMCPt = new TH2F*[fnCuts];
 438                         hESDTruePrimaryMotherInvMassMCPt = new TH2F*[fnCuts];
 439                         hESDTrueSecondaryMotherInvMassPt = new TH2F*[fnCuts];
 440                         hESDTrueSecondaryMotherFromK0sInvMassPt = new TH2F*[fnCuts];
 441                         hESDTrueBckGGInvMassPt = new TH2F*[fnCuts];
 442                         hESDTrueBckContInvMassPt = new TH2F*[fnCuts];
 443                         hESDTrueMotherGGInvMassPt = new TH2F*[fnCuts];
 444                 //}
 445
 446                 for(Int_t iCut = 0; iCut<fnCuts;iCut++){
 447                         TString cutstringElectron =((AliDalitzElectronCuts*)fCutElectronArray->At(iCut))->GetCutNumber();
 448                         TString cutstringMeson= ((AliConversionMesonCuts*)fCutMesonArray->At(iCut))->GetCutNumber();
 449                         TString cutstringGamma = ((AliConversionCuts*)fCutGammaArray->At(iCut))->GetCutNumber();
 450
 451                         fMCList[iCut] = new TList();
 452                         fMCList[iCut]->SetName(Form("%s_%s_%s MC histograms",cutstringGamma.Data(),cutstringElectron.Data(),cutstringMeson.Data()));
 453                         fMCList[iCut]->SetOwner(kTRUE);
 454                         fCutFolder[iCut]->Add(fMCList[iCut]);
 455
 456                                 hMCPi0Pt[iCut] = new TH1F("MC_Pi0_Pt","MC_Pi0_Pt",250,0,25);
 457                                 fMCList[iCut]->Add(hMCPi0Pt[iCut]);
 458                                 hMCEtaPt[iCut] = new TH1F("MC_Eta_Pt","MC_Eta_Pt",250,0,25);
 459                                 fMCList[iCut]->Add(hMCEtaPt[iCut]);
 460                                 hMCPi0InAccPt[iCut] = new TH1F("MC_Pi0DalitzInAcc_Pt","MC_Pi0DalitzInAcc_Pt",250,0,25);
 461                                 fMCList[iCut]->Add(hMCPi0InAccPt[iCut]);
 462                                 hMCEtaInAccPt[iCut] = new TH1F("MC_EtaDalitzInAcc_Pt","MC_EtaDalitzInAcc_Pt",250,0,25);
 463                                 fMCList[iCut]->Add(hMCEtaInAccPt[iCut]);
 464
 465                         fTrueList[iCut] = new TList();
 466                         fTrueList[iCut]->SetName(Form("%s_%s_%s True histograms",cutstringGamma.Data(),cutstringElectron.Data(),cutstringMeson.Data()));
 467                         fTrueList[iCut]->SetOwner(kTRUE);
 468                         fCutFolder[iCut]->Add(fTrueList[iCut]);
 469
 470                         hESDTrueConvGammaPt[iCut] = new TH1F("ESD_TrueConvGamma_Pt","ESD_TrueConvGamma_Pt",250,0,25);
 471                         fTrueList[iCut]->Add(hESDTrueConvGammaPt[iCut]);
 472
 473                         hESDTrueMotherInvMassPt[iCut] = new TH2F("ESD_TrueMother_InvMass_Pt","ESD_TrueMother_InvMass_Pt",1000,0,1,250,0,25);
 474                         fTrueList[iCut]->Add(hESDTrueMotherInvMassPt[iCut]);
 475
 476                         hESDTrueMotherPi0GGInvMassPt[iCut] = new TH2F("ESD_TrueMotherPi0GG_InvMass_Pt","ESD_TrueMotherPi0GG_InvMass_Pt",1000,0,1,250,0,25);
 477                                 fTrueList[iCut]->Add(hESDTrueMotherPi0GGInvMassPt[iCut]);
 478                                 hESDTruePrimaryPi0DalitzESDPtMCPt[iCut] = new TH2F("ESD_TruePrimaryPi0Dalitz_ESDPt_MCPt","ESD_TruePrimaryPi0Dalitz_ESDPt_MCPt",250,0,25,250,0,25);
 479                                 fTrueList[iCut]->Add(hESDTruePrimaryPi0DalitzESDPtMCPt[iCut]);
 480                                 hESDTruePrimaryMotherInvMassMCPt[iCut] = new TH2F("ESD_TruePrimaryMother_InvMass_MCPt","ESD_TrueDalitzPrimaryMother_InvMass_MCPt",1000,0,1,250,0,25);
 481                                 fTrueList[iCut]->Add(hESDTruePrimaryMotherInvMassMCPt[iCut]);
 482                                 hESDTrueSecondaryMotherInvMassPt[iCut] = new TH2F("ESD_TrueDalitzSecondaryMother_InvMass_Pt","ESD_TrueDalitzSecondaryMother_InvMass_Pt",1000,0,1,250,0,25);
 483                                 fTrueList[iCut]->Add(hESDTrueSecondaryMotherInvMassPt[iCut]);
 484 //                              hESDTrueSecondaryMotherFromK0sInvMassPt[iCut] = new TH2F("ESD_TrueDalitzSecondaryMotherFromK0s_InvMass_Pt","ESD_TrueDalitzSecondaryMotherFromK0s_InvMass_Pt",1000,0,1,250,0,25);
 485 //                              fTrueList[iCut]->Add(hESDTrueSecondaryMotherFromK0sInvMassPt[iCut]);
 486                                 hESDTrueBckGGInvMassPt[iCut] = new TH2F("ESD_TrueDalitzBckGG_InvMass_Pt","ESD_TrueDalitzBckGG_InvMass_Pt",1000,0,1,250,0,25);
 487                                 fTrueList[iCut]->Add(hESDTrueBckGGInvMassPt[iCut]);
 488                                 hESDTrueBckContInvMassPt[iCut] = new TH2F("ESD_TrueDalitzBckCont_InvMass_Pt","ESD_TrueDalitzBckCont_InvMass_Pt",1000,0,1,250,0,25);
 489                                 fTrueList[iCut]->Add(hESDTrueBckContInvMassPt[iCut]);
 490 //                              hESDTrueMotherGGInvMassPt[iCut] = new TH2F("ESD_TrueGammaGamma_InvMass_Pt","ESD_TrueGammaGamma_InvMass_Pt",1000,0,1,250,0,25);
 491 //                              fTrueList[iCut]->Add(hESDTrueMotherGGInvMassPt[iCut]);
 492
 493                 }
 494         }
 495
 496         PostData(1, fOutputContainer);
 497
 498 }
 499
 500 //______________________________________________________________________
 501 void AliAnalysisTaskGammaConvDalitzV1::UserExec(Option_t *)
 502 {
 503 //
 504 // Execute analysis for current event
 505 //
 506
 507    fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
 508    if(!fV0Reader){printf("Error: No V0 Reader");return;} // GetV0Reader
 509
 510    Int_t eventQuality = ((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetEventQuality();
 511
 512    if(eventQuality != 0){// Event Not Accepted
 513       for(Int_t iCut = 0; iCut<fnCuts; iCut++){
 514          hNEvents[iCut]->Fill(eventQuality);
 515       }
 516       return;
 517    }
 518
 519    fElecSelector=(AliDalitzElectronSelector*)AliAnalysisManager::GetAnalysisManager()->GetTask("ElectronSelector");
 520    if(!fElecSelector){printf("Error: No ElectronSelector");return;} // GetV0Reader
 521
 522
 523
 524    fMCEvent = MCEvent();
 525    fESDEvent = (AliESDEvent*) InputEvent();
 526    fReaderGammas = fV0Reader->GetReconstructedGammas(); // Gammas from default Cut
 527    fSelectorElectronIndex = fElecSelector->GetReconstructedElectronsIndex(); // Electrons from default Cut
 528    fSelectorPositronIndex = fElecSelector->GetReconstructedPositronsIndex(); // Positrons from default Cut
 529    CountESDTracks(); // Estimate Event Multiplicity
 530
 531    AddTaskContainers(); //Add conatiner
 532
 533     for(Int_t iCut = 0; iCut<fnCuts; iCut++){
 534       fiCut = iCut;
 535
 536       if(fIsHeavyIon && !((AliConversionCuts*)fCutGammaArray->At(iCut))->IsCentralitySelected(fESDEvent)){
 537
 538          hNEvents[iCut]->Fill(1); // Check Centrality --> Not Accepted => eventQuality = 1
 539          continue;
 540
 541       }
 542       hNEvents[iCut]->Fill(eventQuality);
 543
 544       hNGoodESDTracks[iCut]->Fill(fNumberOfESDTracks);
 545
 546       if(fMCEvent){ // Process MC Particle
 547          fMCStack = fMCEvent->Stack();
 548          ProcessMCParticles();
 549       }
 550
 551       ProcessPhotonCandidates(); // Process this cuts gammas
 552       ProcessElectronCandidates(); // Process this cuts gammas
 553       CalculatePi0DalitzCandidates();
 554       CalculateBackground();
 555       UpdateEventByEventData();
 556
 557
 558 //       if(fDoMesonAnalysis){ // Meson Analysis
 559 //          CalculatePi0Candidates(); // Combine Gammas
 560 //          CalculateBackground(); // Combinatorial Background
 561 //          UpdateEventByEventData(); // Store Event for mixed Events
 562 //        }
 563       fGoodGammas->Clear(); // delete this cuts good gammas
 564       fGoodVirtualGammas->Clear(); // delete this cuts good gammas
 565    }
 566
 567     fSelectorElectronIndex.clear();
 568     fSelectorPositronIndex.clear();
 569
 570     PostData( 1, fOutputContainer );
 571 }
 572
 573 void AliAnalysisTaskGammaConvDalitzV1::Terminate(const Option_t *)
 574 {
 575
 576
 577 }
 578
 579 void AliAnalysisTaskGammaConvDalitzV1::AddTaskContainers()
 580 {
 581 //
 582
 583
 584    // cout<<"Entro a dfdsfasdfdsfasdfsdfds terminate: "<<endl;
 585
 586
 587     if( fElecSelector ){
 588
 589          //cout<<"Entro al fElecSelector: "<<endl;
 590
 591          TList *temp0 =  ((AliDalitzElectronCuts*)fElecSelector->GetDalitzElectronCuts())->GetCutHistograms();
 592
 593         if (  temp0 ) {
 594
 595                 if( ! fOutputContainer->Contains( temp0 ) ) {
 596                         fOutputContainer->Add( temp0 );
 597                         cout<<" Add fElecSelector: "<<endl;
 598                 }
 599
 600         }
 601     }
 602
 603
 604     if ( fV0Reader ) {
 605
 606         TList *temp1 = ((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms();
 607
 608         if(  temp1  ) {
 609                 if( ! fOutputContainer->Contains( temp1 ) ) {
 610
 611                         fOutputContainer->Add( temp1 );
 612                         cout<<"Add GammaConv"<<endl;
 613                 }
 614         }
 615     }
 616
 617
 618    // fOutputContainer->Print();
 619
 620    // if( fDebug ) AliInfo("Not to do anything in Terminate");
 621 }
 622
 623
 624 //________________________________________________________________________
 625 void AliAnalysisTaskGammaConvDalitzV1::ProcessPhotonCandidates()
 626 {
 627    Int_t nV0 = 0;
 628    TList *GoodGammasStepOne = new TList();
 629    TList *GoodGammasStepTwo = new TList();
 630    // Loop over Photon Candidates allocated by ReaderV1
 631    for(Int_t i = 0; i < fReaderGammas->GetEntriesFast(); i++){
 632       AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) fReaderGammas->At(i);
 633       if(!PhotonCandidate) continue;
 634       if(!((AliConversionCuts*)fCutGammaArray->At(fiCut))->PhotonIsSelected(PhotonCandidate,fESDEvent)) continue;
 635
 636       if(!((AliConversionCuts*)fCutGammaArray->At(fiCut))->UseElecSharingCut() &&
 637          !((AliConversionCuts*)fCutGammaArray->At(fiCut))->UseToCloseV0sCut()){ // if no post reader loop is required add to events good gammas
 638          fGoodGammas->Add(PhotonCandidate);
 639          hESDConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
 640          if(fMCEvent){
 641             ProcessTruePhotonCandidates(PhotonCandidate);
 642          }
 643       }
 644       else if(((AliConversionCuts*)fCutGammaArray->At(fiCut))->UseElecSharingCut()){ // if Shared Electron cut is enabled, Fill array, add to step one
 645          ((AliConversionCuts*)fCutGammaArray->At(fiCut))->FillElectonLabelArray(PhotonCandidate,nV0);
 646          nV0++;
 647          GoodGammasStepOne->Add(PhotonCandidate);
 648       }
 649       else if(!((AliConversionCuts*)fCutGammaArray->At(fiCut))->UseElecSharingCut() &&
 650               ((AliConversionCuts*)fCutGammaArray->At(fiCut))->UseToCloseV0sCut()){ // shared electron is disabled, step one not needed -> step two
 651          GoodGammasStepTwo->Add(PhotonCandidate);
 652       }
 653    }
 654    if(((AliConversionCuts*)fCutGammaArray->At(fiCut))->UseElecSharingCut()){
 655       for(Int_t i = 0;i<GoodGammasStepOne->GetEntries();i++){
 656          AliAODConversionPhoton *PhotonCandidate= (AliAODConversionPhoton*) GoodGammasStepOne->At(i);
 657          if(!PhotonCandidate) continue;
 658          if(!((AliConversionCuts*)fCutGammaArray->At(fiCut))->RejectSharedElectronV0s(PhotonCandidate,i,GoodGammasStepOne->GetEntries())) continue;
 659          if(!((AliConversionCuts*)fCutGammaArray->At(fiCut))->UseToCloseV0sCut()){ // To Colse v0s cut diabled, step two not needed
 660             fGoodGammas->Add(PhotonCandidate);
 661             hESDConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
 662             if(fMCEvent){
 663                ProcessTruePhotonCandidates(PhotonCandidate);
 664             }
 665          }
 666          else GoodGammasStepTwo->Add(PhotonCandidate); // Close v0s cut enabled -> add to list two
 667       }
 668    }
 669    if(((AliConversionCuts*)fCutGammaArray->At(fiCut))->UseToCloseV0sCut()){
 670       for(Int_t i = 0;i<GoodGammasStepTwo->GetEntries();i++){
 671          AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) GoodGammasStepTwo->At(i);
 672          if(!PhotonCandidate) continue;
 673          if(!((AliConversionCuts*)fCutGammaArray->At(fiCut))->RejectToCloseV0s(PhotonCandidate,GoodGammasStepTwo,i)) continue;
 674          fGoodGammas->Add(PhotonCandidate); // Add gamma to current cut TList
 675          hESDConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt()); // Differences to old V0Reader in p_t due to conversion KF->TLorentzVector
 676          if(fMCEvent){
 677             ProcessTruePhotonCandidates(PhotonCandidate);
 678          }
 679       }
 680    }
 681
 682    delete GoodGammasStepOne;
 683    GoodGammasStepOne = 0x0;
 684    delete GoodGammasStepTwo;
 685    GoodGammasStepTwo = 0x0;
 686 }
 687
 688 //________________________________________________________________________
 689 void AliAnalysisTaskGammaConvDalitzV1::ProcessTruePhotonCandidates(AliAODConversionPhoton *TruePhotonCandidate)
 690 {
 691    // Process True Photons
 692    AliStack *MCStack = fMCEvent->Stack();
 693    TParticle *posDaughter = TruePhotonCandidate->GetPositiveMCDaughter(MCStack);
 694    TParticle *negDaughter = TruePhotonCandidate->GetNegativeMCDaughter(MCStack);
 695
 696    if(posDaughter == NULL || negDaughter == NULL) return; // One particle does not exist
 697    if(posDaughter->GetMother(0) != negDaughter->GetMother(0)){  // Not Same Mother == Combinatorial Bck
 698       return;
 699    }
 700    if(TMath::Abs(posDaughter->GetPdgCode())!=11 || TMath::Abs(negDaughter->GetPdgCode())!=11) return; //One Particle is not electron
 701    if(posDaughter->GetPdgCode()==negDaughter->GetPdgCode()) return; // Same Charge
 702    if(posDaughter->GetUniqueID() != 5 || negDaughter->GetUniqueID() !=5) return;// check if the daughters come from a conversion
 703
 704    TParticle *Photon = TruePhotonCandidate->GetMCParticle(MCStack);
 705    if(Photon->GetPdgCode() != 22) return; // Mother is no Photon
 706
 707    // True Photon
 708    hESDTrueConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
 709 }
 710
 711
 712 //________________________________________________________________________
 713 void AliAnalysisTaskGammaConvDalitzV1::ProcessElectronCandidates(){
 714
 715    Float_t magField = fInputEvent->GetMagneticField();
 716
 717
 718    if( magField  < 0.0 ){
 719        magField =  1.0;
 720    }
 721    else {
 722        magField =  -1.0;
 723    }
 724
 725
 726
 727
 728
 729     vector<Int_t> lGoodElectronIndex(0);
 730     vector<Int_t> lGoodPositronIndex(0);
 731
 732
 733     lGoodElectronIndex.clear();
 734     lGoodPositronIndex.clear();
 735
 736
 737     for(UInt_t i = 0; i < fSelectorElectronIndex.size(); i++){
 738
 739         AliESDtrack* electronCandidate = fESDEvent->GetTrack(fSelectorElectronIndex[i]);
 740         if(! ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->ElectronIsSelected(electronCandidate) ) continue;
 741         lGoodElectronIndex.push_back(   fSelectorElectronIndex[i] );
 742         hESDDalitzElectronPt[fiCut]->Fill(electronCandidate->Pt());
 743
 744     }
 745
 746     for(UInt_t i = 0; i < fSelectorPositronIndex.size(); i++){
 747
 748         AliESDtrack* positronCandidate = fESDEvent->GetTrack(fSelectorPositronIndex[i]);
 749         if(! ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->ElectronIsSelected(positronCandidate) ) continue;
 750         lGoodPositronIndex.push_back(   fSelectorPositronIndex[i] );
 751         hESDDalitzPositronPt[fiCut]->Fill(positronCandidate->Pt());
 752
 753     }
 754
 755
 756
 757
 758
 759
 760    for(UInt_t i = 0; i < lGoodElectronIndex.size(); i++){
 761
 762
 763                 AliESDtrack* electronCandidate = fESDEvent->GetTrack(lGoodElectronIndex[i]);
 764
 765                 AliKFParticle electronCandidateKF( *electronCandidate->GetConstrainedParam(), ::kElectron );
 766
 767                 for(UInt_t j = 0; j < lGoodPositronIndex.size(); j++){
 768                     AliESDtrack* positronCandidate = fESDEvent->GetTrack(lGoodPositronIndex[j]);
 769                          AliKFParticle positronCandidateKF( *positronCandidate->GetConstrainedParam(), ::kPositron );
 770
 771                         if( ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->DoPsiPairCut() == kTRUE ){
 772
 773
 774
 775
 776                             Double_t psiPair = ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->GetPsiPair(electronCandidate,positronCandidate);
 777                             Double_t deltaPhi = magField * TVector2::Phi_mpi_pi( electronCandidate->GetConstrainedParam()->Phi()-positronCandidate->GetConstrainedParam()->Phi());
 778
 779
 780
 781                             if( ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->IsFromGammaConversion(psiPair,deltaPhi ) ){
 782                                  continue;
 783                             }
 784                               hESDEposEnegPsiPairDPhi[fiCut]->Fill(deltaPhi,psiPair);
 785
 786                           }
 787
 788
 789                           AliKFConversionPhoton* virtualPhoton = new AliKFConversionPhoton(electronCandidateKF,positronCandidateKF);
 790
 791                          //if(fUseImprovedVertex == kTRUE){
 792                                 AliKFVertex primaryVertexImproved(*fInputEvent->GetPrimaryVertex());
 793                                 primaryVertexImproved+=*virtualPhoton;
 794                                 virtualPhoton->SetProductionVertex(primaryVertexImproved);
 795                          //}
 796
 797                           virtualPhoton->SetTrackLabels( lGoodPositronIndex[j], lGoodElectronIndex[i]);
 798
 799                           if(fMCEvent){
 800
 801                                // AliStack *fMCStack= fMCEvent->Stack();
 802                                 Int_t labeln=TMath::Abs(electronCandidate->GetLabel());
 803                                 Int_t labelp=TMath::Abs(positronCandidate->GetLabel());
 804                                 TParticle *fNegativeMCParticle = fMCStack->Particle(labeln);
 805                                 TParticle *fPositiveMCParticle = fMCStack->Particle(labelp);
 806                                 if( fPositiveMCParticle && fNegativeMCParticle) {
 807                                     virtualPhoton->SetMCLabelPositive(labelp);
 808                                     virtualPhoton->SetMCLabelNegative(labeln);
 809                                 }
 810                           }
 811                           AliAODConversionPhoton *vphoton = new AliAODConversionPhoton(virtualPhoton); //To Apply PsiPairCut
 812                           fGoodVirtualGammas->Add(  vphoton );
 813             }
 814    }
 815
 816 }
 817
 818 //________________________________________________________________________
 819 void AliAnalysisTaskGammaConvDalitzV1::CalculatePi0DalitzCandidates(){
 820
 821    // Conversion Gammas
 822
 823
 824
 825
 826    if( fGoodGammas->GetEntries() > 0 && fGoodVirtualGammas->GetEntries() > 0 ){
 827
 828
 829       for(Int_t GammaIndex=0; GammaIndex<fGoodGammas->GetEntries(); GammaIndex++){
 830
 831          AliAODConversionPhoton *gamma=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(GammaIndex));
 832
 833          for(Int_t virtualGammaIndex=0;virtualGammaIndex<fGoodVirtualGammas->GetEntries();virtualGammaIndex++){
 834
 835             AliAODConversionPhoton *Vgamma=dynamic_cast<AliAODConversionPhoton*>(fGoodVirtualGammas->At(virtualGammaIndex));
 836             //Check for same Electron ID
 837             if(gamma->GetTrackLabelPositive() == Vgamma->GetTrackLabelPositive() ||
 838                gamma->GetTrackLabelNegative() == Vgamma->GetTrackLabelNegative() ||
 839                gamma->GetTrackLabelNegative() == Vgamma->GetTrackLabelPositive() ||
 840                gamma->GetTrackLabelPositive() == Vgamma->GetTrackLabelNegative() ) continue;
 841
 842             AliAODConversionMother *pi0cand = new AliAODConversionMother(gamma,Vgamma);
 843             pi0cand->SetLabels(GammaIndex,virtualGammaIndex);
 844
 845               if( ( ((AliConversionMesonCuts*)fCutMesonArray->At(fiCut))->MesonIsSelected(pi0cand,kTRUE)) ){
 846                hESDMotherInvMassPt[fiCut]->Fill(pi0cand->M(),pi0cand->Pt());
 847               // if(pi0cand->GetAlpha()<0.1){
 848               //    hESDMotherInvMassEalpha[fiCut]->Fill(pi0cand->M(),pi0cand->E());
 849               // }
 850                Int_t zbin= fBGHandler[fiCut]->GetZBinIndex(fESDEvent->GetPrimaryVertex()->GetZ());
 851                Int_t mbin = 0;
 852                if(((AliConversionCuts*)fCutGammaArray->At(fiCut))->UseTrackMultiplicity()){
 853                   mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fNumberOfESDTracks);
 854                } else {
 855                   mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGoodGammas->GetEntries());
 856                }
 857                Double_t sparesFill[4] = {pi0cand->M(),pi0cand->Pt(),zbin,mbin};
 858                sESDMotherInvMassPtZM[fiCut]->Fill(sparesFill,1);
 859                if(fMCEvent){
 860                   ProcessTrueMesonCandidates(pi0cand,gamma,Vgamma);
 861                }
 862               }
 863             delete pi0cand;
 864             pi0cand=0x0;
 865          }
 866       }
 867     }
 868 }
 869
 870 //________________________________________________________________________
 871 void AliAnalysisTaskGammaConvDalitzV1::CalculateBackground(){
 872
 873    Int_t zbin= fBGHandler[fiCut]->GetZBinIndex(fESDEvent->GetPrimaryVertex()->GetZ());
 874    Int_t mbin = 0;
 875
 876    Int_t method = 0;
 877
 878    method = ((AliDalitzElectronCuts*) fCutElectronArray->At(fiCut))->GetBKGMethod();
 879
 880
 881    if(((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->UseTrackMultiplicity()){
 882       mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fNumberOfESDTracks);
 883    } else {
 884       mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGoodGammas->GetEntries());
 885    }
 886
 887    if( method == 1 || method == 2 ) {
 888
 889       AliGammaConversionAODBGHandler::GammaConversionVertex *bgEventVertex = NULL;
 890
 891       if( ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->UseTrackMultiplicity() ) {
 892
 893          for(Int_t nEventsInBG=0;nEventsInBG<fBGHandler[fiCut]->GetNBGEvents();nEventsInBG++){
 894
 895             AliGammaConversionAODVector *previousEventV0s = fBGHandler[fiCut]->GetBGGoodV0s(zbin,mbin,nEventsInBG);
 896
 897             if(fMoveParticleAccordingToVertex == kTRUE && method == 1){
 898                bgEventVertex = fBGHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG);
 899             }
 900
 901             for(Int_t iCurrent=0;iCurrent<fGoodVirtualGammas->GetEntries();iCurrent++){
 902                AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGoodVirtualGammas->At(iCurrent));
 903
 904                 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
 905                   AliAODConversionPhoton previousGoodV0 = (AliAODConversionPhoton)(*(previousEventV0s->at(iPrevious)));
 906
 907                   if(fMoveParticleAccordingToVertex == kTRUE && method == 1 ){
 908                      MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
 909                   }
 910
 911                   AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&previousGoodV0);
 912                   if( ( ((AliConversionMesonCuts*)fCutMesonArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE) ) ){
 913                      hESDMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
 914                      Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),zbin,mbin};
 915                      sESDMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
 916                   }
 917                   delete backgroundCandidate;
 918                   backgroundCandidate = 0x0;
 919                }
 920             }
 921          }
 922       }
 923       else{
 924          for(Int_t nEventsInBG=0;nEventsInBG <fBGHandler[fiCut]->GetNBGEvents();nEventsInBG++){
 925             AliGammaConversionAODVector *previousEventV0s = fBGHandler[fiCut]->GetBGGoodV0s(zbin,mbin,nEventsInBG);
 926             if(previousEventV0s){
 927                if(fMoveParticleAccordingToVertex == kTRUE && method == 1){
 928                   bgEventVertex = fBGHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG);
 929                }
 930                for(Int_t iCurrent=0;iCurrent<fGoodVirtualGammas->GetEntries();iCurrent++){
 931                   AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGoodVirtualGammas->At(iCurrent));
 932                   for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
 933
 934                      AliAODConversionPhoton previousGoodV0 = (AliAODConversionPhoton)(*(previousEventV0s->at(iPrevious)));
 935
 936                      if(fMoveParticleAccordingToVertex == kTRUE && method ==1){
 937                         MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
 938                      }
 939
 940                      AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&previousGoodV0);
 941
 942                      if((((AliConversionMesonCuts*)fCutMesonArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE))){
 943                         hESDMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
 944                         Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),zbin,mbin};
 945                         sESDMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
 946                      }
 947                      delete backgroundCandidate;
 948                      backgroundCandidate = 0x0;
 949                   }
 950                }
 951             }
 952          }
 953       }
 954    }
 955
 956    else if( method == 3 ){
 957
 958             for(Int_t iCurrent=0;iCurrent<fGoodVirtualGammas->GetEntries();iCurrent++){
 959
 960                AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGoodVirtualGammas->At(iCurrent));
 961
 962                 for(Int_t iPrevious=0;iPrevious<fGammasPool[fiCut]->GetEntries();iPrevious++){
 963
 964                   AliAODConversionPhoton previousGoodV0 = *(AliAODConversionPhoton*)((fGammasPool[fiCut]->At(iPrevious) ));
 965
 966
 967                   AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&previousGoodV0);
 968
 969                   if((((AliConversionMesonCuts*)fCutMesonArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE))){
 970
 971                      hESDMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
 972                      Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),zbin,mbin};
 973                      sESDMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
 974                   }
 975                   delete backgroundCandidate;
 976                   backgroundCandidate = 0x0;
 977                }
 978             }
 979    }
 980
 981 }
 982 //________________________________________________________________________
 983 void AliAnalysisTaskGammaConvDalitzV1::UpdateEventByEventData(){
 984    //see header file for documentation
 985
 986    Int_t method = 0;
 987
 988    method = ( (AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->GetBKGMethod();
 989
 990    if( method == 1 ) {
 991
 992    if(fGoodGammas->GetEntries() >0 ){
 993
 994       if(((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->UseTrackMultiplicity()){
 995          fBGHandler[fiCut]->AddEvent(fGoodGammas,fESDEvent->GetPrimaryVertex()->GetX(),fESDEvent->GetPrimaryVertex()->GetY(),fESDEvent->GetPrimaryVertex()->GetZ(),fNumberOfESDTracks);
 996       }
 997
 998       else{ // means we use #V0s for multiplicity
 999          fBGHandler[fiCut]->AddEvent(fGoodGammas,fESDEvent->GetPrimaryVertex()->GetX(),fESDEvent->GetPrimaryVertex()->GetY(),fESDEvent->GetPrimaryVertex()->GetZ(),fGoodGammas->GetEntries());
1000       }
1001     }
1002    }
1003
1004    else if ( method == 2 ){
1005
1006     if(fGoodVirtualGammas->GetEntries() >0 ){
1007       if(((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->UseTrackMultiplicity()){
1008          fBGHandler[fiCut]->AddEvent(fGoodVirtualGammas,fESDEvent->GetPrimaryVertex()->GetX(),fESDEvent->GetPrimaryVertex()->GetY(),fESDEvent->GetPrimaryVertex()->GetZ(),fNumberOfESDTracks);
1009       }
1010       else{ // means we use #V0s for multiplicity
1011          fBGHandler[fiCut]->AddEvent(fGoodVirtualGammas,fESDEvent->GetPrimaryVertex()->GetX(),fESDEvent->GetPrimaryVertex()->GetY(),fESDEvent->GetPrimaryVertex()->GetZ(),fGoodGammas->GetEntries());
1012      }
1013    }
1014   }
1015   else if ( method  == 3 ) {
1016
1017                for(Int_t index = 0; index < fGoodGammas->GetEntries(); index ++){
1018
1019                             if ( fGammasPool[fiCut]->GetEntries() > 10 ){
1020                                     fGammasPool[fiCut]->RemoveLast();
1021                             }
1022                  fGammasPool[fiCut]->AddFirst(dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(index)));
1023                }
1024       }
1025 }
1026 //______________________________________________________________________
1027 void AliAnalysisTaskGammaConvDalitzV1::ProcessTrueMesonCandidates(AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate, AliAODConversionPhoton *TrueVirtualGammaCandidate)
1028 {
1029
1030    // Process True Mesons
1031
1032    AliStack *MCStack = fMCEvent->Stack();
1033
1034
1035    if(TrueGammaCandidate->GetV0Index()<fESDEvent->GetNumberOfV0s()){
1036
1037
1038       Bool_t isTruePi0 = kFALSE;
1039       Bool_t isTrueEta = kFALSE;
1040       Int_t gammaMCLabel = TrueGammaCandidate->GetMCParticleLabel(MCStack);
1041       Int_t gammaMotherLabel = -1;
1042
1043
1044
1045
1046       if(gammaMCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1047
1048
1049          // Daughters Gamma 0
1050          TParticle * negativeMC = (TParticle*)TrueGammaCandidate->GetNegativeMCDaughter(MCStack);
1051          TParticle * positiveMC = (TParticle*)TrueGammaCandidate->GetPositiveMCDaughter(MCStack);
1052          TParticle * gammaMC = (TParticle*)MCStack->Particle(gammaMCLabel);
1053
1054
1055          if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1056
1057             if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
1058
1059                     if(gammaMC->GetPdgCode() == 22){ // ... with Gamma Mother
1060                      gammaMotherLabel=gammaMC->GetFirstMother();
1061                     }
1062            }
1063         }
1064        }
1065
1066
1067          Int_t virtualGammaMCLabel = TrueVirtualGammaCandidate->GetMCParticleLabel(MCStack);
1068          Int_t virtualGammaMotherLabel = -1;
1069          Int_t virtualGamma = 1;
1070
1071          if(virtualGammaMCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1072             // Daughters Gamma 1
1073             TParticle * negativeMC = (TParticle*)TrueVirtualGammaCandidate->GetNegativeMCDaughter(MCStack);
1074             TParticle * positiveMC = (TParticle*)TrueVirtualGammaCandidate->GetPositiveMCDaughter(MCStack);
1075             TParticle * virtualGammaMotherMC = (TParticle*)MCStack->Particle(virtualGammaMCLabel);
1076
1077             if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1078
1079             if( virtualGammaMotherMC->GetPdgCode() != 22 ){
1080                   virtualGammaMotherLabel=virtualGammaMCLabel;
1081             }
1082
1083              else if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
1084
1085                     virtualGammaMotherLabel=virtualGammaMotherMC->GetFirstMother();
1086                     virtualGamma = 0; //no virtual gamma
1087
1088              }
1089             }
1090          }
1091
1092
1093          if(gammaMotherLabel>=0 && ( gammaMotherLabel == virtualGammaMotherLabel) ){
1094
1095             if(((TParticle*)MCStack->Particle(virtualGammaMotherLabel))->GetPdgCode() == 111){
1096                isTruePi0=kTRUE;
1097             }
1098
1099             if(((TParticle*)MCStack->Particle(virtualGammaMotherLabel))->GetPdgCode() == 221){
1100                isTrueEta=kTRUE;
1101             }
1102
1103          }
1104
1105          if(isTruePi0 || isTrueEta ){ // True Pion or Eta
1106             if ( virtualGamma == 1 ) { //True Dalitz
1107             hESDTrueMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1108             if(gammaMotherLabel <= MCStack->GetNprimary()){ // Only primary pi0 for efficiency calculation
1109                hESDTruePrimaryMotherInvMassMCPt[fiCut]->Fill(Pi0Candidate->M(),((TParticle*)MCStack->Particle(virtualGammaMotherLabel))->Pt());
1110                if(isTruePi0){ // Only primaries for unfolding
1111                   hESDTruePrimaryPi0DalitzESDPtMCPt[fiCut]->Fill(Pi0Candidate->Pt(),((TParticle*)MCStack->Particle(virtualGammaMotherLabel))->Pt());
1112                }
1113             }
1114             if(gammaMotherLabel > MCStack->GetNprimary()){ // Secondary Meson
1115                hESDTrueSecondaryMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1116                //if (((TParticle*)MCStack->Particle(virtualGammaMotherLabel))->GetMother(0) >-1){
1117                //   if(MCStack->Particle(((TParticle*)MCStack->Particle(virtualGammaMotherLabel))->GetMother(0))->GetPdgCode()==kK0Short){
1118                //      hESDTrueSecondaryMotherFromK0sInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1119                //   }
1120                }
1121             }
1122             else if ( virtualGamma == 0 ){
1123               hESDTrueMotherPi0GGInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt()); // Pi0 from GG
1124             }
1125          }
1126
1127         if(!isTruePi0 && !isTrueEta){ // Background
1128             if(gammaMotherLabel>-1 && virtualGammaMotherLabel>-1 && virtualGamma == 0){ // Both Tracks are Photons and have a mother but not Pi0 or Eta
1129                hESDTrueBckGGInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1130             } else { // No photon or without mother
1131                hESDTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1132             }
1133         }
1134     }
1135 }
1136 //________________________________________________________________________
1137 void AliAnalysisTaskGammaConvDalitzV1::MoveParticleAccordingToVertex(AliAODConversionPhoton* particle,const AliGammaConversionAODBGHandler::GammaConversionVertex *vertex){
1138    //see header file for documentation
1139
1140    Double_t dx = vertex->fX - fESDEvent->GetPrimaryVertex()->GetX();
1141    Double_t dy = vertex->fY - fESDEvent->GetPrimaryVertex()->GetY();
1142    Double_t dz = vertex->fZ - fESDEvent->GetPrimaryVertex()->GetZ();
1143
1144    Double_t movedPlace[3] = {particle->GetConversionX() - dx,particle->GetConversionY() - dy,particle->GetConversionZ() - dz};
1145    particle->SetConversionPoint(movedPlace);
1146 }
1147
1148
1149 //________________________________________________________________________
1150 void AliAnalysisTaskGammaConvDalitzV1::CountESDTracks(){
1151
1152    AliESDtrackCuts *EsdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
1153    // Using standard function for setting Cuts
1154    Bool_t selectPrimaries=kTRUE;
1155    EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
1156    EsdTrackCuts->SetMaxDCAToVertexZ(2);
1157    EsdTrackCuts->SetEtaRange(-0.8, 0.8);
1158    EsdTrackCuts->SetPtRange(0.15);
1159
1160    fNumberOfESDTracks = 0;
1161    for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
1162       AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
1163       if(!curTrack) continue;
1164       if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
1165    }
1166    delete EsdTrackCuts;
1167    EsdTrackCuts=0x0;
1168
1169    return;
1170 }
1171
1172 //_____________________________________________________________________________
1173 void AliAnalysisTaskGammaConvDalitzV1::ProcessMCParticles()
1174 {
1175
1176
1177    // Loop over all primary MC particle
1178    for(Int_t i = 0; i < fMCStack->GetNprimary(); i++) {
1179       TParticle* particle = (TParticle *)fMCStack->Particle(i);
1180       if (!particle) continue;
1181
1182
1183          if(((AliConversionMesonCuts*)fCutMesonArray->At(fiCut))->MesonIsSelectedMCDalitz(particle,fMCStack)){
1184             if(particle->GetPdgCode() == 111)hMCPi0Pt[fiCut]->Fill(particle->Pt()); // All MC Pi0
1185             if(particle->GetPdgCode() == 221)hMCEtaPt[fiCut]->Fill(particle->Pt()); // All MC Eta
1186             // Check the acceptance for both gammas
1187             if(particle->GetNDaughters() == 3){
1188                TParticle* gamma    = 0;
1189                TParticle* electron = 0;
1190                TParticle* positron = 0;
1191
1192
1193                for(Int_t index=particle->GetFirstDaughter(); index<= particle->GetLastDaughter();index++){
1194
1195
1196                     TParticle* temp = (TParticle*)fMCStack->Particle( index );
1197
1198                     switch( temp->GetPdgCode() ) {
1199                                 case ::kPositron:
1200                                         electron = temp;
1201                                         break;
1202                                 case ::kElectron:
1203                                         positron = temp;
1204                                         break;
1205                                 case ::kGamma:
1206                                         gamma    = temp;
1207                                         break;
1208                         }
1209                }
1210                if( gamma  && electron && positron ) {
1211
1212                  if(((AliConversionCuts*)fCutGammaArray->At(fiCut))->PhotonIsSelectedMC(gamma,fMCStack,kFALSE) &&
1213                     TMath::Abs( electron->Eta() ) < ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->GetEtaCut()  &&
1214                     TMath::Abs( positron->Eta() ) < ((AliDalitzElectronCuts*)fCutElectronArray->At(fiCut))->GetEtaCut() ){
1215                   if(particle->GetPdgCode() == 111)hMCPi0InAccPt[fiCut]->Fill(particle->Pt()); // MC Pi0Dalitz with gamma and e+e- in acc
1216                   if(particle->GetPdgCode() == 221)hMCEtaInAccPt[fiCut]->Fill(particle->Pt()); // MC EtaDalitz with gamma and e+e- in acc
1217                }
1218               }
1219             }
1220          }
1221       }
1222 }