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