PWG4/GammaConv/AliAnalysisTaskGammaConvDalitz.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 "AliGammaConversionHistograms.h"
  47 #include "AliV0Reader.h"
  48 #include "AliKFVertex.h"
  49
  50 #include "AliAnalysisTaskGammaConvDalitz.h"
  51 #include "TH1.h"
  52
  53 ClassImp( AliAnalysisTaskGammaConvDalitz )
  54
  55 //-----------------------------------------------------------------------------------------------
  56 AliAnalysisTaskGammaConvDalitz::AliAnalysisTaskGammaConvDalitz():
  57  AliAnalysisTaskSE(),
  58  fStack(0),
  59  fGCMCEvent(0),
  60  fESDEvent(0),
  61  fEposCandidateIndex(),
  62  fEnegCandidateIndex(),
  63  fGammaCandidatePosIndex(),
  64  fGammaCandidateNegIndex(),
  65  fGammaCandidates(0),
  66  fGammaPool(0),
  67  fPoolMaxSize(10),
  68  fGamPoolPos(0),
  69  fBGEventHandler(0),
  70  fOutputContainer(0),
  71  fMCTruth(0),
  72  fV0Reader(0),
  73  fESDpid(0),
  74  fESDtrackCuts(0),
  75  fITSsaTrackCuts(0),
  76  fESDpidCuts(0),
  77  fHistograms(0),
  78  fStandalone(kFALSE),
  79  fDoMC(kFALSE),
  80  fComputeBkg(kTRUE),
  81  fUseBayesPID(kFALSE),
  82  fUseTrackIndexCut(kTRUE),
  83  fUsePsiPairCut(kTRUE),
  84  fUseMassCut(kFALSE),
  85  fUseGammaCut(kFALSE),
  86  fReadMagFieldSign(kTRUE),
  87  fUseAliKF(kFALSE),
  88  fMagFieldSign(1),
  89  fkElectronMass(0.00051099891),
  90  fPsiPairCut(0.45),
  91  fDeltaPhiCutMin(0.),
  92  fDeltaPhiCutMax(0.12),
  93  fMassCutMin(0.),
  94  fMassCutMax(0.1),
  95  fNSigmaBelowElecTPCbethe(-2.),
  96  fNSigmaAboveElecTPCbethe(3.),
  97  fNSigmaAbovePionTPCbethe(3.),
  98  fNSigmaAboveKaonTPCbethe(3.),
  99  fNSigmaAboveProtonTPCbethe(3.),
 100  fTrkSelectionCriteria(kGlobalTrack)
 101 {
 102 //
 103 // Default constructor
 104 //
 105         AdoptITSsaTrackCuts();
 106         AdoptESDtrackCuts();
 107         AdoptESDpidCuts();
 108
 109         fGammaPool = new TClonesArray("AliKFParticle", fPoolMaxSize);
 110         fGammaPool->SetOwner(kTRUE);
 111 }
 112
 113 //-----------------------------------------------------------------------------------------------
 114 AliAnalysisTaskGammaConvDalitz::AliAnalysisTaskGammaConvDalitz( const char* name ):
 115  AliAnalysisTaskSE( name ),
 116  fStack(0),
 117  fGCMCEvent(0),
 118  fESDEvent(0),
 119  fEposCandidateIndex(),
 120  fEnegCandidateIndex(),
 121  fGammaCandidatePosIndex(),
 122  fGammaCandidateNegIndex(),
 123  fGammaCandidates(0),
 124  fGammaPool(0),
 125  fPoolMaxSize(10),
 126  fGamPoolPos(0),
 127  fBGEventHandler(0),
 128  fOutputContainer(0),
 129  fMCTruth(0),
 130  fV0Reader(0),
 131  fESDpid(0),
 132  fESDtrackCuts(0),
 133  fITSsaTrackCuts(0),
 134  fESDpidCuts(0),
 135  fHistograms(0),
 136  fStandalone(kFALSE),
 137  fDoMC(kFALSE),
 138  fComputeBkg(kTRUE),
 139  fUseBayesPID(kFALSE),
 140  fUseTrackIndexCut(kTRUE),
 141  fUsePsiPairCut(kTRUE),
 142  fUseMassCut(kFALSE),
 143  fUseGammaCut(kFALSE),
 144  fReadMagFieldSign(kTRUE),
 145  fUseAliKF(kFALSE),
 146  fMagFieldSign(1),
 147  fkElectronMass(0.00051099891),
 148  fPsiPairCut(0.45),
 149  fDeltaPhiCutMin(0.),
 150  fDeltaPhiCutMax(0.12),
 151  fMassCutMin(0.),
 152  fMassCutMax(0.1),
 153  fNSigmaBelowElecTPCbethe(-2.),
 154  fNSigmaAboveElecTPCbethe(3.),
 155  fNSigmaAbovePionTPCbethe(3.),
 156  fNSigmaAboveKaonTPCbethe(3.),
 157  fNSigmaAboveProtonTPCbethe(3.),
 158  fTrkSelectionCriteria(kGlobalTrack)
 159 {
 160     // Common I/O in slot 0
 161     DefineInput (0, TChain::Class());
 162
 163     // Your private output
 164     DefineOutput(1, TList::Class());
 165 //  DefineOutput(2, AliCFContainer::Class());  // for CF
 166
 167     AdoptITSsaTrackCuts();
 168     AdoptESDtrackCuts();
 169     AdoptESDpidCuts();
 170    // fkElectronMass = TDatabasePDG::Instance()->GetParticle( ::kElectron )->Mass(); //
 171
 172     fGammaPool = new TClonesArray("AliKFParticle", fPoolMaxSize);
 173     fGammaPool->SetOwner(kTRUE);
 174 }
 175
 176 //-----------------------------------------------------------------------------------------------
 177 AliAnalysisTaskGammaConvDalitz::~AliAnalysisTaskGammaConvDalitz()
 178 {
 179 //
 180 // virtual destructor
 181 //
 182
 183     if( fOutputContainer )          delete fOutputContainer;
 184     if( fHistograms )               delete fHistograms;
 185     if( fStandalone && fV0Reader )  delete fV0Reader;
 186     if( fITSsaTrackCuts )           delete fITSsaTrackCuts;
 187     if( fESDtrackCuts )             delete fESDtrackCuts;
 188     if( fESDpidCuts )               delete fESDpidCuts;
 189     if( fGammaCandidates)           delete fGammaCandidates;
 190     if( fGammaPool )                delete fGammaPool;
 191 }
 192
 193 //-----------------------------------------------------------------------------------------------
 194 void AliAnalysisTaskGammaConvDalitz::ConnectInputData(Option_t *option)
 195 {
 196 //
 197 // Connect Input Data
 198 //
 199     if( fDebug ) AliInfo("=> ConnectInputData");
 200
 201     AliAnalysisTaskSE::ConnectInputData(option);
 202
 203     if( fV0Reader == 0 )
 204     {
 205         AliFatal("There is not pointer to AliV0Reader object!!!");
 206     }
 207 }
 208
 209 //-----------------------------------------------------------------------------------------------
 210 void AliAnalysisTaskGammaConvDalitz::UserCreateOutputObjects()
 211 {
 212 //
 213 // Create ouput objects
 214 //
 215     if( fDebug ) AliInfo("=> UserCreateOutputObjects");
 216
 217     // Create the output container
 218     if( fOutputContainer != 0 )
 219     {
 220         delete fOutputContainer;
 221     }
 222
 223     fOutputContainer = new TList();
 224
 225     // Add the histograms to the output container
 226     fHistograms->GetOutputContainer( fOutputContainer );
 227     fOutputContainer->SetOwner(kTRUE);
 228
 229     PostData( 1, fOutputContainer );
 230 }
 231
 232 //-----------------------------------------------------------------------------------------------
 233 void AliAnalysisTaskGammaConvDalitz::UserExec(Option_t */*option*/)
 234 {
 235 //
 236 // Execute analysis for current event
 237 //
 238     if( fDebug ) AliInfo("=> UserExec");
 239
 240     if( fV0Reader == 0 )
 241     {
 242         AliFatal("no pointer to AliV0Reader");
 243         return;
 244     }
 245
 246     // Create list of gamma candidates in standalone mode
 247     // otherwise use the created ones by AliAnalysisTaskGammaConversion
 248     if( fStandalone )
 249     {
 250
 251        AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 252        AliESDInputHandler *esdHandler=0;
 253         if ( (esdHandler=dynamic_cast<AliESDInputHandler*>(man->GetInputEventHandler())) && esdHandler->GetESDpid() ){
 254             AliV0Reader::SetESDpid(esdHandler->GetESDpid());
 255      } else {
 256     //load esd pid bethe bloch parameters depending on the existance of the MC handler
 257     // yes: MC parameters
 258     // no:  data parameters
 259         if (!AliV0Reader::GetESDpid()){
 260           if (MCEvent() ) {
 261                 AliV0Reader::InitESDpid();
 262             } else {
 263                 AliV0Reader::InitESDpid(1);
 264             }
 265         }
 266       }
 267
 268
 269         if (MCEvent() ) {
 270
 271     // To avoid crashes due to unzip errors. Sometimes the trees are not there.
 272         AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 273
 274         if (!mcHandler){
 275          AliError("Could not retrive MC event handler!");
 276          return;
 277         }
 278
 279         if (!mcHandler->InitOk() ){
 280
 281         return;
 282         }
 283         if (!mcHandler->TreeK() ){
 284         return;
 285         }
 286          if (!mcHandler->TreeTR() ) {
 287         return;
 288          }
 289       }
 290
 291
 292
 293
 294        fV0Reader->SetInputAndMCEvent( InputEvent(), MCEvent() );
 295        fV0Reader->Initialize();
 296     }
 297
 298     if( fV0Reader->CheckForPrimaryVertex() == kFALSE )
 299     {
 300         if( fDebug ) AliInfo("no contributors to primary vertex");
 301         return;
 302     }
 303
 304     if( fV0Reader->CheckForPrimaryVertexZ() == kFALSE  )
 305     {
 306
 307         if( fDebug ) AliInfo("z vertex out of range");
 308         return;
 309     }
 310
 311     // Get Pointers
 312    fBGEventHandler = fV0Reader->GetBGHandler();
 313    fESDpid = fV0Reader->GetESDpid();
 314    fESDEvent = fV0Reader->GetESDEvent();
 315    if(fDoMC && MCEvent())
 316    {
 317         fStack= MCEvent()->Stack();
 318         fGCMCEvent=MCEvent();
 319    }
 320
 321     // Read the magnetic field sign from ESD
 322     if ( fReadMagFieldSign == kTRUE )
 323     {
 324         fMagFieldSign = (fESDEvent->GetMagneticField() < 0) ? 1 : -1;
 325     }
 326
 327     // Process MC information
 328     if(fDoMC)
 329     {
 330         ProcessMCData();
 331     }
 332
 333     if( fStandalone ){
 334             while(fV0Reader->NextV0()){}; //SelectGammas
 335             fV0Reader->ResetV0IndexNumber();
 336     }
 337
 338     CreateListOfDalitzPairCandidates();
 339     ProcessGammaElectronsForDalitzAnalysis();
 340
 341     if ( fStandalone ){
 342
 343       fV0Reader->UpdateEventByEventData();
 344
 345     }
 346
 347     PostData( 1, fOutputContainer );
 348 }
 349
 350
 351 void AliAnalysisTaskGammaConvDalitz::Terminate(Option_t */*option*/)
 352 {
 353 //
 354     if( fDebug ) AliInfo("Not to do anything in Terminate");
 355 }
 356
 357 //-----------------------------------------------------------------------------------------------
 358 void AliAnalysisTaskGammaConvDalitz::AdoptITSsaTrackCuts( AliESDtrackCuts* esdCuts )
 359 {
 360 //
 361 // set user ITSsa track cuts
 362 //
 363     if( fITSsaTrackCuts ) delete fITSsaTrackCuts;
 364
 365     if( esdCuts )
 366     {
 367         fITSsaTrackCuts = esdCuts;
 368     }
 369     else
 370     {
 371         // default cuts
 372         fITSsaTrackCuts = new AliESDtrackCuts("Default ITSsa track cuts for Pi0 Dalitz decay");
 373
 374         fITSsaTrackCuts->SetEtaRange( -0.9, 0.9 );
 375         fITSsaTrackCuts->SetAcceptKinkDaughters(kFALSE);
 376
 377         fITSsaTrackCuts->SetMinNClustersITS(2);
 378         fITSsaTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
 379         fITSsaTrackCuts->SetRequireITSRefit(kTRUE);
 380
 381         fITSsaTrackCuts->SetRequireSigmaToVertex(kTRUE);
 382         fITSsaTrackCuts->SetMaxNsigmaToVertex(3);
 383
 384         fITSsaTrackCuts->SetRequireITSStandAlone(kTRUE);
 385     }
 386 }
 387
 388 //-----------------------------------------------------------------------------------------------
 389 void AliAnalysisTaskGammaConvDalitz::AdoptESDtrackCuts( AliESDtrackCuts* esdCuts )
 390 {
 391 //
 392 // set user global track cuts
 393 //
 394     if( fESDtrackCuts ) delete fESDtrackCuts;
 395
 396     if( esdCuts )
 397     {
 398         fESDtrackCuts = esdCuts;
 399     }
 400     else
 401     {
 402         //default cuts
 403         fESDtrackCuts = new AliESDtrackCuts("Default global track cuts for Pi0 Dalitz decay");
 404
 405         fESDtrackCuts->SetEtaRange( -0.9, 0.9 );
 406         fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
 407
 408         fESDtrackCuts->SetMinNClustersITS(2);
 409         fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
 410         fESDtrackCuts->SetRequireITSRefit(kTRUE);
 411
 412         fESDtrackCuts->SetRequireSigmaToVertex(kTRUE);
 413         fESDtrackCuts->SetMaxNsigmaToVertex(3);
 414
 415         fESDtrackCuts->SetMinNClustersTPC(80);
 416         fESDtrackCuts->SetMaxChi2PerClusterTPC(4.);
 417         fESDtrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
 418         fESDtrackCuts->SetRequireTPCRefit(kTRUE);
 419     }
 420 }
 421
 422 //-----------------------------------------------------------------------------------------------
 423 void AliAnalysisTaskGammaConvDalitz::AdoptESDpidCuts( AliESDpidCuts* esdPIDCuts )
 424 {
 425 //
 426 // set user pid cuts
 427 //
 428     if( fESDpidCuts ) delete fESDpidCuts;
 429     if( esdPIDCuts )
 430     {
 431         fESDpidCuts = esdPIDCuts;
 432     }
 433     else // default cuts
 434     {
 435         fESDpidCuts = new AliESDpidCuts("Electrons", "Electron PID cuts");
 436      //   fESDpidCuts->SetTPCnSigmaCut(AliPID::kElectron, 3.);
 437         fESDpidCuts->SetTPCnSigmaCut(AliPID::kElectron, -4., 6.);
 438     }
 439 }
 440
 441 //-----------------------------------------------------------------------------------------------
 442 void AliAnalysisTaskGammaConvDalitz::ProcessMCData()
 443 {
 444 //
 445 // Process generation
 446 //
 447         if( fDebug ) AliInfo("=> ProcessMCData");
 448
 449         fHistograms->FillTable("Table_Generation", 0);  //number of events
 450
 451         for ( Int_t i = 0; i < fStack->GetNtrack(); i++ )
 452         {
 453                 TParticle* iParticle = fStack->Particle( i );
 454                 if( !iParticle ) continue;
 455
 456                 if ( i >= fStack->GetNprimary() ) continue;  // is primary?
 457                 if ( iParticle->GetPdgCode() != ::kPi0 ) continue;  // is Pi0?
 458
 459                 if( iParticle->GetNDaughters() == 2 &&
 460                     fStack->Particle(iParticle->GetFirstDaughter())->GetPdgCode() == ::kGamma &&
 461                     fStack->Particle(iParticle->GetLastDaughter())->GetPdgCode() == ::kGamma )
 462                 {
 463                         fHistograms->FillTable("Table_Generation", 1);  // pi0 -> gg
 464                 }
 465
 466                 if ( iParticle->GetNDaughters() != 3 ) continue;    // Num == 3 (e+,e-,gamma)
 467
 468                 // Check for Pi0 Dalitz decay
 469                 TParticle* eposPi0 = 0;
 470                 TParticle* enegPi0 = 0;
 471                 TParticle* gammaPi0 = 0;
 472
 473                 for( Int_t idxPi0 = iParticle->GetFirstDaughter(); idxPi0 <= iParticle->GetLastDaughter(); idxPi0++ )
 474                 {
 475                         switch(fStack->Particle(idxPi0)->GetPdgCode())
 476                         {
 477                                 case ::kPositron:
 478                                         eposPi0 = fStack->Particle(idxPi0);
 479                                         break;
 480                                 case ::kElectron:
 481                                         enegPi0 = fStack->Particle(idxPi0);
 482                                         break;
 483                                 case ::kGamma:
 484                                         gammaPi0 = fStack->Particle(idxPi0);
 485                                         break;
 486                         }
 487                 }
 488
 489                 if (eposPi0==0 || enegPi0==0 || gammaPi0==0) continue;
 490
 491                 // found a Pi0 Dalitz decay
 492
 493                 fHistograms->FillTable("Table_Generation", 2);
 494                 fHistograms->FillHistogram("MC_Pi0Dalitz_P", iParticle->P());
 495                 fHistograms->FillHistogram("MC_Pi0Dalitz_Pt", iParticle->Pt());
 496                 fHistograms->FillHistogram("MC_Pi0Dalitz_Eta", iParticle->Eta());
 497                 fHistograms->FillHistogram("MC_Pi0Dalitz_Pt_vs_Y", Rapidity(iParticle), iParticle->Pt());
 498                 fHistograms->FillHistogram("MC_EposDalitz_Pt", eposPi0->Pt());
 499                 fHistograms->FillHistogram("MC_EposDalitz_Eta", eposPi0->Eta());
 500                 fHistograms->FillHistogram("MC_EnegDalitz_Pt", enegPi0->Pt());
 501                 fHistograms->FillHistogram("MC_EnegDalitz_Eta", enegPi0->Eta());
 502                 fHistograms->FillHistogram("MC_GammaPi0Dalitz_Pt", gammaPi0->Pt());
 503                 fHistograms->FillHistogram("MC_GammaPi0Dalitz_Eta", gammaPi0->Eta());
 504
 505                 // Angle between the gamma and the plane e+e-
 506                 TVector3 ePosMom( eposPi0->Px(), eposPi0->Py(), eposPi0->Pz() );
 507                 TVector3 eNegMom( enegPi0->Px(), enegPi0->Py(), enegPi0->Pz() );
 508                 TVector3 gamMom( gammaPi0->Px(), gammaPi0->Py() , gammaPi0->Pz() );
 509                 TVector3 planeEposEneg =  eNegMom.Cross( ePosMom );
 510                 Double_t anglePlaneGamma = planeEposEneg.Angle(gamMom);
 511
 512                 fHistograms->FillHistogram("MC_EposEnegDalitz_Angle", ePosMom.Angle(eNegMom) );
 513
 514                 fHistograms->FillHistogram("MC_EposEnegDalitz_GammaPi0_Angle", anglePlaneGamma);
 515                 fHistograms->FillHistogram("MC_EposEnegDalitz_GammaPi0_Angle_vs_P", anglePlaneGamma, gammaPi0->P());
 516                 fHistograms->FillHistogram("MC_EposEnegDalitz_GammaPi0_Angle_vs_Pt", anglePlaneGamma, gammaPi0->Pt());
 517
 518
 519         // check for gamma conversion
 520         Bool_t daugGammaElectron    = kFALSE;
 521         Bool_t daugGammaPositron    = kFALSE;  // acceptance
 522         Bool_t daugGammaElectronAll = kFALSE;
 523         Bool_t daugGammaPositronAll = kFALSE;
 524
 525         // is the gamma converted? -> has 2 daughter e+e-
 526         // are e+ e- from gamma in the acceptance for the V0s
 527         if( gammaPi0->GetNDaughters() >= 2 )
 528         {
 529             for( Int_t tIndex=gammaPi0->GetFirstDaughter(); tIndex<=gammaPi0->GetLastDaughter(); ++tIndex )
 530             {
 531                 TParticle* tmpDaughter = fStack->Particle(tIndex);
 532
 533                 if( tmpDaughter->GetUniqueID() != kPPair ) continue; // check if the daughters come from a conversion
 534                 if( tmpDaughter->GetPdgCode() == ::kElectron )
 535                 { // e+
 536                     daugGammaElectronAll = kTRUE;
 537
 538                     if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() &&
 539                         ((TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()) < tmpDaughter->R() &&
 540                         (tmpDaughter->R()< fV0Reader->GetMaxRCut() ) )
 541                     {
 542                         daugGammaElectron = kTRUE;
 543                     }
 544                 }
 545                 else if( tmpDaughter->GetPdgCode() == ::kPositron )
 546                 {
 547                      daugGammaPositronAll = kTRUE;
 548                     if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() &&
 549                         ((TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()) < tmpDaughter->R() &&
 550                         (tmpDaughter->R() < fV0Reader->GetMaxRCut() ) )
 551                     {
 552                         daugGammaPositron = kTRUE;
 553                     }
 554                 }
 555             }
 556         }
 557
 558
 559        if(  daugGammaElectronAll && daugGammaPositronAll )
 560        {
 561           TParticle* tmpDaughter = fStack->Particle( gammaPi0->GetFirstDaughter() );
 562           fHistograms->FillHistogram("MC_GC_GammaPi0Dalitz_All_Z_vs_R",tmpDaughter->Vz(),tmpDaughter->R() );
 563        }
 564
 565         Float_t  etaMin, etaMax;
 566         fESDtrackCuts->GetEtaRange( etaMin, etaMax );
 567
 568         // e+e- pair within acceptance
 569         if ( TMath::Abs( eposPi0->Eta() ) < etaMax  && TMath::Abs( enegPi0->Eta() ) < etaMax )
 570         {
 571             fHistograms->FillHistogram("MC_Acceptance_EposDalitz_Pt", eposPi0->Pt());
 572             fHistograms->FillHistogram("MC_Acceptance_EposDalitz_Eta", eposPi0->Eta());
 573             fHistograms->FillHistogram("MC_Acceptance_EnegDalitz_Pt", enegPi0->Pt());
 574             fHistograms->FillHistogram("MC_Acceptance_EnegDalitz_Eta", enegPi0->Eta());
 575             fHistograms->FillHistogram("MC_Acceptance_DalitzPair_EposPt_vs_EnegPt", eposPi0->Pt(), enegPi0->Pt());
 576         }
 577
 578         // Pi0 (e+e-gamma) within acceptance
 579         //cout<<"Gamma Eta Cut"<<fV0Reader->GetEtaCut()<<endl;
 580
 581         if ( ( TMath::Abs( gammaPi0->Eta() ) < fV0Reader->GetEtaCut() && gammaPi0->R() < fV0Reader->GetMaxRCut() ) &&
 582              TMath::Abs( eposPi0->Eta() ) < etaMax  && TMath::Abs( enegPi0->Eta() ) < etaMax )
 583         {
 584             fHistograms->FillTable("Table_Generation",3);  //
 585             fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Pt",iParticle->Pt());
 586             fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Eta",iParticle->Eta());
 587             fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Radius",iParticle->R());
 588             fHistograms->FillHistogram("MC_Acceptance_GammaPi0Dalitz_Pt",gammaPi0->Pt());
 589             fHistograms->FillHistogram("MC_Acceptance_GammaPi0Dalitz_Eta",gammaPi0->Eta());
 590             fHistograms->FillHistogram("MC_Acceptance_EposPi0Dalitz_Pt",eposPi0->Pt());
 591             fHistograms->FillHistogram("MC_Acceptance_EposPi0Dalitz_Eta",eposPi0->Eta());
 592             fHistograms->FillHistogram("MC_Acceptance_EnegPi0Dalitz_Pt",enegPi0->Pt());
 593             fHistograms->FillHistogram("MC_Acceptance_EnegPi0Dalitz_Eta",enegPi0->Eta());
 594             fHistograms->FillHistogram("MC_Acceptance_DalitzPair_OpeningAngle", ePosMom.Angle(eNegMom) );
 595             fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Pt_vs_Y", Rapidity(iParticle), iParticle->Pt());
 596
 597            // Pi0 within acceptance with gamma converted
 598
 599             if ( daugGammaElectron && daugGammaPositron )
 600             {
 601                 fHistograms->FillTable("Table_Generation",4); //
 602
 603                 fHistograms->FillHistogram("MC_Acceptance_GC_Pi0Dalitz_Pt",iParticle->Pt());
 604                 fHistograms->FillHistogram("MC_Acceptance_GC_Pi0Dalitz_Eta",iParticle->Eta());
 605                 fHistograms->FillHistogram("MC_Acceptance_GC_EposPi0Dalitz_Pt",eposPi0->Pt());
 606                 fHistograms->FillHistogram("MC_Acceptance_GC_EposPi0Dalitz_Eta",eposPi0->Eta());
 607                 fHistograms->FillHistogram("MC_Acceptance_GC_EnegPi0Dalitz_Pt",enegPi0->Pt());
 608                 fHistograms->FillHistogram("MC_Acceptance_GC_EnegPi0Dalitz_Eta",enegPi0->Eta());
 609                 fHistograms->FillHistogram("MC_Acceptance_GC_GammaPi0Dalitz_Pt",gammaPi0->Pt());
 610                 fHistograms->FillHistogram("MC_Acceptance_GC_GammaPi0Dalitz_Eta",gammaPi0->Eta());
 611                 //fHistograms->FillHistogram("MC_Acceptance_GC_Gamma_Angle",anglePlaneGamma);
 612                 //fHistograms->FillHistogram("MC_Acceptance_GC_Gamma_Angle_vs_Pt",anglePlaneGamma,gammaPi0->Pt());
 613                 TParticle* tmpDaughter = fStack->Particle( gammaPi0->GetFirstDaughter() );
 614                 fHistograms->FillHistogram("MC_Acceptance_GC_GammaPi0Dalitz_Z_vs_R",tmpDaughter->Vz(),tmpDaughter->R() );
 615                 fHistograms->FillHistogram("MC_Acceptance_GC_Pi0Dalitz_Pt_vs_Y", Rapidity(iParticle), iParticle->Pt());
 616             }
 617         }
 618     }
 619 }
 620
 621 //-----------------------------------------------------------------------------------------------
 622 void AliAnalysisTaskGammaConvDalitz::CreateListOfDalitzPairCandidates()
 623 {
 624 //
 625 // Dalitz pair candidates
 626 //
 627         if( fDebug ) AliInfo("=> CreateListOfDalitzPairCandidates");
 628
 629         fEposCandidateIndex.clear();
 630         fEnegCandidateIndex.clear();
 631
 632         fHistograms->FillTable("Table_Cuts", 0);
 633
 634         for( Int_t i = 0; i < fESDEvent->GetNumberOfTracks(); ++i )
 635         {
 636                 AliESDtrack* iTrack = fESDEvent->GetTrack(i);
 637                 if ( !iTrack ) continue;
 638
 639
 640                 Double_t p[3];
 641
 642                 if ( !iTrack->GetConstrainedPxPyPz(p) ) continue;
 643
 644                 TVector3 iMom(p[0],p[1],p[2]);
 645
 646                 //
 647                 // Check track cuts and find track type
 648                 //
 649
 650                 Bool_t isTrackAccepted = 0;
 651                 Int_t trackType = -1;
 652                 switch(fTrkSelectionCriteria)
 653                 {
 654                         case kITSsaTrack:
 655                                 isTrackAccepted = fITSsaTrackCuts->AcceptTrack( iTrack );
 656                                 trackType = kITSsaTrack;
 657                                 break;
 658
 659                         case kGlobalTrack:
 660                                 isTrackAccepted = fESDtrackCuts->AcceptTrack( iTrack );
 661                                 trackType = kGlobalTrack;
 662                                 break;
 663
 664                         case kITSsaGlobalTrack:
 665                                 if(fITSsaTrackCuts->AcceptTrack( iTrack ) || fESDtrackCuts->AcceptTrack( iTrack ))
 666                                 {
 667                                         isTrackAccepted = kTRUE;
 668                                         if(fITSsaTrackCuts->AcceptTrack( iTrack )) trackType = kITSsaTrack;
 669                                         else trackType = kGlobalTrack;
 670                                 }
 671                                 break;
 672                 }
 673
 674                 if(!isTrackAccepted) continue;
 675
 676                 //
 677                 // PID
 678                 //
 679
 680                 Int_t pid=-1;
 681                 Int_t pidMC=-1;
 682
 683                 if(fUseBayesPID)
 684                 {
 685                         pid = GetBayesPid(iTrack,trackType);
 686                 }
 687                 else
 688                 {
 689                         pid = GetNSigmaPid(iTrack,trackType);
 690                 }
 691
 692                 if( fDoMC )
 693                 {
 694                         pidMC = GetMonteCarloPid(iTrack);
 695                         // pid table
 696                         Int_t iLabel = TMath::Abs(iTrack->GetLabel());
 697                         TParticle* iParticle = fStack->Particle(iLabel);
 698                         FillPidTable(iParticle, pid);
 699                 }
 700
 701                 // ITS standalone tracks
 702                 if( trackType == kITSsaTrack)
 703                 {
 704                         Double_t mom = iTrack->GetP();
 705                         Double_t signal = iTrack->GetITSsignal();
 706
 707                         fHistograms->FillHistogram( "ESD_ITSsa_dEdx_vs_P", mom, signal );
 708
 709                         if( pid == AliPID::kElectron )
 710                         {
 711
 712                                 fHistograms->FillHistogram( "ESD_ITSsa_PidCut_dEdx_vs_P", mom, signal );
 713                                 if(fDoMC && pid == pidMC)
 714                                 {
 715                                         fHistograms->FillHistogram( "MC_ESD_ITSsa_PidCut_dEdx_vs_P", mom, signal );
 716                                 }
 717                         }
 718
 719                         if( fDoMC && pidMC == AliPID::kElectron)
 720                         {
 721                                 fHistograms->FillHistogram( "MC_ESD_ITSsa_Electron_dEdx_vs_P", mom, signal );
 722                         }
 723                 }
 724
 725                 else  // global tracks
 726                 {
 727                         const AliExternalTrackParam *in = iTrack->GetInnerParam();
 728                         Double_t mom = in->GetP();
 729                         Double_t signal = iTrack->GetTPCsignal();
 730
 731                         fHistograms->FillHistogram( "ESD_TPC_dEdx_vs_P", mom, signal );
 732
 733                         if( fDoMC && pidMC == AliPID::kElectron )
 734                         {
 735                                 fHistograms->FillHistogram( "MC_ESD_TPC_Electron_dEdx_vs_P", mom, signal );
 736                         }
 737
 738                         if( pid == AliPID::kElectron )
 739                         {
 740                                 fHistograms->FillHistogram( "ESD_TPC_PidCut_dEdx_vs_P", mom, signal );
 741                                 if(fDoMC && pid == pidMC)
 742                                 {
 743                                         fHistograms->FillHistogram( "MC_ESD_TPC_PidCut_dEdx_vs_P", mom, signal );
 744                                 }
 745                         }
 746                 }
 747
 748                 if( AliPID::kElectron != pid) continue;
 749
 750                 // electron track candidates from here
 751
 752                 if( iTrack->GetSign() > 0 )
 753                 {
 754                         fEposCandidateIndex.push_back(i);
 755                 }
 756                 else
 757                 {
 758                         fEnegCandidateIndex.push_back(i);
 759                 }
 760         }
 761
 762         // gamma candidates
 763         GetGammaCandidates(fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
 764
 765         if(fDoMC)
 766         {
 767                 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
 768                 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(1.,(Double_t)pi0Dalitz->GetEntriesFast());
 769                 delete pi0Dalitz;
 770         }
 771
 772         if(fUseTrackIndexCut) // remove repeated tracks
 773         {
 774                 ESDtrackIndexCut(fEposCandidateIndex,fEnegCandidateIndex, fGammaCandidates);
 775
 776                 if(fDoMC)
 777                 {
 778                         TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
 779                         ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(2.,(Double_t)pi0Dalitz->GetEntriesFast());
 780                         delete pi0Dalitz;
 781                 }
 782         }
 783
 784         if(fUsePsiPairCut) // remove electrons from gamma conversions
 785         {
 786                 PsiPairCut(fEposCandidateIndex,fEnegCandidateIndex);
 787
 788                 if(fDoMC)
 789                 {
 790                         TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
 791                         ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(3.,(Double_t)pi0Dalitz->GetEntriesFast());
 792                         delete pi0Dalitz;
 793                 }
 794         }
 795
 796         if( fUseMassCut )
 797         {
 798                 MassCut(fEposCandidateIndex, fEnegCandidateIndex);
 799
 800                 if(fDoMC)
 801                 {
 802                         TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
 803                         ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(4.,(Double_t)pi0Dalitz->GetEntriesFast());
 804                         delete pi0Dalitz;
 805                 }
 806         }
 807
 808         if(fUseGammaCut)
 809         {
 810                 AngleEposEnegGammaCut(fEposCandidateIndex,fEnegCandidateIndex,fV0Reader->GetCurrentEventGoodV0s(), fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
 811
 812                 if(fDoMC)
 813                 {
 814                         TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
 815                         ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(5.,(Double_t)pi0Dalitz->GetEntriesFast());
 816                         delete pi0Dalitz;
 817                 }
 818         }
 819 }
 820
 821 //-----------------------------------------------------------------------------------------------
 822 void AliAnalysisTaskGammaConvDalitz::ProcessGammaElectronsForDalitzAnalysis()
 823 {
 824 //
 825 // Process gamma and electrons for pi0 Dalitz decay
 826 //
 827         if( fDebug ) AliInfo("=> ProcessGammaElectronsForDalitzAnalysis");
 828
 829         fHistograms->FillTable( "Table_Reconstruction", 0); // number of events
 830
 831         TClonesArray* ePosCandidates = IndexToAliKFParticle(fEposCandidateIndex, ::kPositron);
 832
 833         for(Int_t i=0; i < ePosCandidates->GetEntriesFast(); ++i)
 834         {
 835                 AliKFParticle* epos = (AliKFParticle*) ePosCandidates->At(i);
 836                 fHistograms->FillHistogram("ESD_EposCandidates_Pt", epos->GetPt());
 837                 fHistograms->FillHistogram("ESD_EposCandidates_Eta", epos->GetEta());
 838                 fHistograms->FillTable( "Table_Reconstruction", 1);
 839         }
 840
 841         TClonesArray* eNegCandidates = IndexToAliKFParticle(fEnegCandidateIndex, ::kElectron);
 842
 843         for(Int_t i=0; i < eNegCandidates->GetEntriesFast(); ++i)
 844         {
 845                 AliKFParticle* eneg = (AliKFParticle*) eNegCandidates->At(i);
 846                 fHistograms->FillHistogram("ESD_EnegCandidates_Pt", eneg->GetPt());
 847                 fHistograms->FillHistogram("ESD_EnegCandidates_Eta", eneg->GetEta());
 848                 fHistograms->FillTable( "Table_Reconstruction", 2);
 849         }
 850
 851         TClonesArray* dalitzPairCandidates = FindDalitzPair(ePosCandidates, eNegCandidates);
 852         for(Int_t i=0; i < dalitzPairCandidates->GetEntriesFast(); ++i)
 853         {
 854                 TLorentzVector* dalitz = (TLorentzVector*)dalitzPairCandidates->At(i);
 855
 856                 fHistograms->FillHistogram("ESD_DalitzPairCandidates_Pt", dalitz->Pt());
 857                 fHistograms->FillHistogram("ESD_DalitzPairCandidates_Mass", dalitz->M());
 858         }
 859
 860         // gamma candidates
 861         for(Int_t i=0; i < fGammaCandidates->GetEntriesFast(); ++i)
 862         {
 863                 AliKFParticle* gamma = (AliKFParticle*) fGammaCandidates->At(i);
 864                 fHistograms->FillHistogram("ESD_GammaCandidates_Pt", gamma->GetPt());
 865                 fHistograms->FillHistogram("ESD_GammaCandidates_Eta", gamma->GetEta());
 866         }
 867
 868         // psi pair for all candidates
 869         //if(fUsePsiPairCut)
 870         FillPsiPair(ePosCandidates,eNegCandidates,"ESD_EposEneg_PsiPair_vs_DPhi");
 871
 872         // Angle epos,eneg gamma
 873         FillAngle(ePosCandidates, fGammaCandidates, "ESD_EposEneg_GammaCandidates_Angle");
 874         FillAngle(eNegCandidates, fGammaCandidates, "ESD_EposEneg_GammaCandidates_Angle");
 875
 876         TClonesArray* pi0Candidates = FindPi0Dalitz(ePosCandidates, eNegCandidates, fGammaCandidates);
 877         for(Int_t i=0; i < pi0Candidates->GetEntriesFast(); ++i)
 878         {
 879                 TLorentzVector* pi0 = (TLorentzVector*)pi0Candidates->At(i);
 880
 881                 fHistograms->FillHistogram("ESD_Pi0_P", pi0->P());
 882                 fHistograms->FillHistogram("ESD_Pi0_Pt", pi0->Pt());
 883                 fHistograms->FillHistogram("ESD_Pi0_Eta", pi0->Eta());
 884                 fHistograms->FillHistogram("ESD_Pi0_Y", pi0->Rapidity());
 885                 fHistograms->FillHistogram("ESD_Pi0_Phi", pi0->Phi());
 886                 fHistograms->FillHistogram("ESD_Pi0_Pt_vs_Y", pi0->Rapidity(), pi0->Pt());
 887                 fHistograms->FillHistogram("ESD_Pi0_Mass", pi0->M());
 888                 fHistograms->FillHistogram("ESD_Pi0_Mass_vs_Pt", pi0->Pt(), pi0->M());
 889                 fHistograms->FillHistogram("ESD_Pi0_Mass_vs_Y",  pi0->Rapidity(),pi0->M());
 890                 fHistograms->FillHistogram("ESD_Pi0_Mass_vs_Eta", pi0->Eta(), pi0->M());
 891         }
 892
 893         delete dalitzPairCandidates;
 894         delete pi0Candidates;
 895
 896         if(fComputeBkg)
 897         {
 898                 // 1) e+e- with with gammas used in the signal
 899                 TClonesArray* pi0Bkg = FindPi0Dalitz(ePosCandidates, eNegCandidates, fGammaPool);
 900
 901                 for(Int_t i=0; i < pi0Bkg->GetEntriesFast(); ++i)
 902                 {
 903                         TLorentzVector* pi0 = (TLorentzVector*)pi0Bkg->At(i);
 904                         fHistograms->FillHistogram("ESD_BKG_PrevGamma", pi0->M());
 905                         fHistograms->FillHistogram("ESD_BKG_PrevGamma_vs_Pt", pi0->Pt(), pi0->M());
 906                 }
 907
 908                 if(ePosCandidates->GetEntriesFast() > 0 &&
 909                    eNegCandidates->GetEntriesFast() > 0 &&
 910                    fGammaCandidates->GetEntriesFast() > 0)
 911                 {
 912                         UpdateGammaPool(fGammaCandidates);
 913                 }
 914
 915                 delete pi0Bkg;
 916
 917                 // 2) e+e- with gammas from a pool of events
 918                 TClonesArray* gammaBGHandler = GammasFromBGHandler();
 919                 pi0Bkg = FindPi0Dalitz(ePosCandidates, eNegCandidates, gammaBGHandler);
 920
 921                 for(Int_t i=0; i < pi0Bkg->GetEntriesFast(); ++i)
 922                 {
 923                         TLorentzVector* pi0 = (TLorentzVector*)pi0Bkg->At(i);
 924                         fHistograms->FillHistogram("ESD_BKG_BGHandler", pi0->M());
 925                         fHistograms->FillHistogram("ESD_BKG_BGHandler_vs_Pt", pi0->Pt(), pi0->M());
 926                 }
 927
 928                 delete pi0Bkg;
 929
 930                 // 3) e+ with e-, gamma from a pool of events
 931                 TClonesArray* elecBGHandler = ElectronFromBGHandler();
 932                 pi0Bkg = FindPi0Dalitz(ePosCandidates, elecBGHandler, gammaBGHandler);
 933
 934                 for(Int_t i=0; i < pi0Bkg->GetEntriesFast(); ++i)
 935                 {
 936                         TLorentzVector* pi0 = (TLorentzVector*)pi0Bkg->At(i);
 937                         fHistograms->FillHistogram("ESD_BKG_Electron", pi0->M());
 938                         fHistograms->FillHistogram("ESD_BKG_Electron_vs_Pt", pi0->Pt(), pi0->M());
 939                 }
 940
 941                 if(eNegCandidates->GetEntriesFast() > 0)
 942                 {
 943                         UpdateElectronPool(eNegCandidates);
 944                 }
 945
 946                 delete gammaBGHandler;
 947                 delete elecBGHandler;
 948                 delete pi0Bkg;
 949
 950         }
 951
 952         delete ePosCandidates;
 953         delete eNegCandidates;
 954
 955         if(fDoMC)
 956         {
 957                 TClonesArray* ePosPi0Dalitz = FindElectronFromPi0Dalitz(fEposCandidateIndex, ::kPositron);
 958                 for(Int_t i=0; i < ePosPi0Dalitz->GetEntriesFast(); ++i)
 959                 {
 960                         AliKFParticle* epos = (AliKFParticle*) ePosPi0Dalitz->At(i);
 961                         fHistograms->FillHistogram("MC_ESD_EposDalitz_Pt", epos->GetPt());
 962                         fHistograms->FillHistogram("MC_ESD_EposDalitz_Eta", epos->GetEta());
 963                         fHistograms->FillTable( "Table_Reconstruction", 3);
 964                 }
 965
 966                 TClonesArray* eNegPi0Dalitz = FindElectronFromPi0Dalitz(fEnegCandidateIndex, ::kElectron);
 967                 for(Int_t i=0; i < eNegPi0Dalitz->GetEntriesFast(); ++i)
 968                 {
 969                         AliKFParticle* eneg = (AliKFParticle*) eNegPi0Dalitz->At(i);
 970                         fHistograms->FillHistogram("MC_ESD_EnegDalitz_Pt", eneg->GetPt());
 971                         fHistograms->FillHistogram("MC_ESD_EnegDalitz_Eta", eneg->GetEta());
 972                         fHistograms->FillTable( "Table_Reconstruction", 4);
 973                 }
 974
 975                 TClonesArray* dalitzPair = FindDalitzPair(fEposCandidateIndex, fEnegCandidateIndex);
 976                 for(Int_t i=0; i < dalitzPair->GetEntriesFast(); ++i)
 977                 {
 978                         TLorentzVector* dalitz = (TLorentzVector*) dalitzPair->At(i);
 979                         fHistograms->FillHistogram("MC_ESD_DalitzPair_Pt", dalitz->Pt());
 980                         fHistograms->FillHistogram("MC_ESD_DalitzPair_Mass", dalitz->M());
 981                         fHistograms->FillHistogram( "Table_Reconstruction", 5 );
 982                 }
 983
 984                 // psi pair for dalitz pairs
 985                 //if(fUsePsiPairCut)
 986                 FillPsiPair(ePosPi0Dalitz,eNegPi0Dalitz,"MC_ESD_DalitzPair_PsiPair_vs_DPhi");
 987
 988                 delete ePosPi0Dalitz;
 989                 delete eNegPi0Dalitz;
 990                 delete dalitzPair;
 991
 992                 // all gammas
 993                 TClonesArray* gamma = FindGamma(fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
 994                 for(Int_t i=0; i < gamma->GetEntriesFast(); ++i)
 995                 {
 996                         AliKFParticle* iGamma = (AliKFParticle*) gamma->At(i);
 997                         fHistograms->FillHistogram("MC_ESD_Gamma_Pt", iGamma->GetPt());
 998                         fHistograms->FillHistogram("MC_ESD_Gamma_Eta", iGamma->GetEta());
 999                 }
1000
1001                 delete gamma;
1002
1003                 // gamma from pi0 dalitz
1004                 TClonesArray* gammaPi0Dalitz = FindGammaFromPi0Dalitz(fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
1005                 for(Int_t i=0; i < gammaPi0Dalitz->GetEntriesFast(); ++i)
1006                 {
1007                         AliKFParticle* iGamma = (AliKFParticle*) gammaPi0Dalitz->At(i);
1008                         fHistograms->FillHistogram("MC_ESD_GammaPi0Dalitz_Pt", iGamma->GetPt());
1009                         fHistograms->FillHistogram("MC_ESD_GammaPi0Dalitz_Eta", iGamma->GetEta());
1010                         fHistograms->FillTable( "Table_Reconstruction", 6);
1011                 }
1012
1013                 delete gammaPi0Dalitz;
1014
1015                 TClonesArray* pi0Dalitz = FindPi0Dalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
1016
1017                 for(Int_t i=0; i < pi0Dalitz->GetEntriesFast(); ++i)
1018                 {
1019                         TLorentzVector* pi0 = (TLorentzVector*) pi0Dalitz->At(i);
1020
1021                         fHistograms->FillHistogram("MC_ESD_Pi0_P", pi0->P());
1022                         fHistograms->FillHistogram("MC_ESD_Pi0_Pt", pi0->Pt());
1023                         fHistograms->FillHistogram("MC_ESD_Pi0_Eta", pi0->Eta());
1024                         fHistograms->FillHistogram("MC_ESD_Pi0_Y", pi0->Rapidity());
1025                         fHistograms->FillHistogram("MC_ESD_Pi0_Phi", pi0->Phi());
1026                         fHistograms->FillHistogram("MC_ESD_Pi0_Pt_vs_Y", pi0->Rapidity(), pi0->Pt());
1027                         fHistograms->FillHistogram("MC_ESD_Pi0_Mass", pi0->M());
1028                         fHistograms->FillHistogram("MC_ESD_Pi0_Mass_vs_Pt", pi0->Pt(), pi0->M());
1029                         fHistograms->FillHistogram("MC_ESD_Pi0_Mass_vs_Y",  pi0->Rapidity(),pi0->M());
1030                         fHistograms->FillHistogram("MC_ESD_Pi0_Mass_vs_Eta", pi0->Eta(), pi0->M());
1031                         fHistograms->FillHistogram( "Table_Reconstruction", 7);
1032                 }
1033                 delete pi0Dalitz;
1034
1035                 // psi pair for electrons from gamma conversions assuming they came from main vertex
1036                 // if(fUsePsiPairCut)
1037                 for(UInt_t i=0; i < fEposCandidateIndex.size(); ++i)
1038                 {
1039                         AliESDtrack* posTrack = fESDEvent->GetTrack(fEposCandidateIndex[i]);
1040                         Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1041
1042                         for(UInt_t j=0; j < fEnegCandidateIndex.size(); ++j)
1043                         {
1044                                 AliESDtrack* negTrack = fESDEvent->GetTrack(fEnegCandidateIndex[j]);
1045                                 Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1046
1047                                 if(!IsFromGammaConversion(posLabel,negLabel)) continue;
1048
1049                                 Double_t psiPair = GetPsiPair(posTrack, negTrack);
1050                                 Double_t deltaPhi = fMagFieldSign * TVector2::Phi_mpi_pi( negTrack->GetConstrainedParam()->Phi()-posTrack->GetConstrainedParam()->Phi());
1051
1052                                 fHistograms->FillHistogram("MC_ESD_EposEnegGamma_PsiPair_vs_DPhi", deltaPhi, psiPair);
1053                         }
1054                 }
1055                 // FIXME: eta -> e+e-gamma
1056         }
1057 }
1058
1059 //--------------------------------------------------------------------------
1060 Double_t AliAnalysisTaskGammaConvDalitz::Rapidity(const TParticle* p) const
1061 {
1062 //
1063 // Get rapidity
1064 //
1065         const double kEPSILON=1.e-16;
1066
1067         if(p->Energy() - TMath::Abs(p->Pz()) < kEPSILON )
1068         {
1069                 return 1.e10;
1070         }
1071         return 0.5*TMath::Log( (p->Energy()+p->Pz()) / (p->Energy()-p->Pz()) );
1072 }
1073
1074 //--------------------------------------------------------------------------
1075 void AliAnalysisTaskGammaConvDalitz::FillPsiPair(const TClonesArray* pos, const TClonesArray* neg, const TString& hName)
1076 {
1077 //
1078 // Fill histogram with psipair(pos,neg)
1079 //
1080         for(Int_t i=0; i < pos->GetEntriesFast(); ++i )
1081         {
1082                 AliKFParticle* posKF = (AliKFParticle*) pos->At(i);
1083                 for( Int_t j=0; j < neg->GetEntriesFast(); ++j )
1084                 {
1085                         AliKFParticle* negKF = (AliKFParticle*) neg->At(j);
1086                         Double_t psiPair = GetPsiPair(posKF, negKF);
1087                         Double_t deltaPhi = fMagFieldSign * TVector2::Phi_mpi_pi( negKF->GetPhi() - posKF->GetPhi());
1088                         fHistograms->FillHistogram(hName, deltaPhi, psiPair);
1089                 }
1090         }
1091 }
1092
1093 //--------------------------------------------------------------------------
1094 void AliAnalysisTaskGammaConvDalitz::FillAngle(const TClonesArray* x, const TClonesArray* y, const TString& hName)
1095 {
1096 //
1097 // Fill histogram with angle(x,y)
1098 //
1099         for(Int_t i=0; i < x->GetEntriesFast(); ++i )
1100         {
1101                 AliKFParticle* xKF = (AliKFParticle*) x->At(i);
1102                 TVector3 xMom(xKF->Px(),xKF->Py(),xKF->Pz());
1103                 for( Int_t j=0; j < y->GetEntriesFast(); ++j )
1104                 {
1105                         AliKFParticle* yKF = (AliKFParticle*) y->At(j);
1106                         TVector3 yMom(yKF->Px(),yKF->Py(),yKF->Pz());
1107                         fHistograms->FillHistogram(hName, xMom.Angle(yMom));
1108                 }
1109         }
1110 }
1111
1112 //--------------------------------------------------------------------------
1113 void AliAnalysisTaskGammaConvDalitz::FillPidTable(const TParticle* p, Int_t pid)
1114 {
1115 //
1116 // Fill table with pid info
1117 //
1118         Int_t iGen=-1;
1119         switch(TMath::Abs(p->GetPdgCode()))
1120         {
1121                 case ::kElectron:   iGen=0; break;
1122                 case ::kMuonMinus:  iGen=1; break;
1123                 case ::kPiPlus:     iGen=2; break;
1124                 case ::kKPlus:      iGen=3; break;
1125                 case ::kProton:     iGen=4; break;
1126                 default: iGen=-1;
1127         }
1128
1129         int jRec=-1;
1130         if(pid > -1 && pid < 5) jRec = pid;
1131
1132         if ((iGen > -1) && (jRec > -1))
1133         {
1134                 fHistograms->FillTable("Table_PID", iGen, jRec);
1135                 // sum
1136                 fHistograms->FillTable("Table_PID", iGen, 5);
1137                 fHistograms->FillTable("Table_PID", 5, jRec);
1138         }
1139 }
1140
1141 //--------------------------------------------------------------------------
1142 void AliAnalysisTaskGammaConvDalitz::GetGammaCandidates(TClonesArray*& gamma, vector<Int_t>& posIndex, vector<Int_t>& negIndex)
1143 {
1144 //
1145 // Make a copy of gamma candidates from V0reader
1146 //
1147         posIndex.clear();
1148         negIndex.clear();
1149
1150         if(gamma) delete gamma;
1151
1152         TClonesArray* gammaV0 = fV0Reader->GetCurrentEventGoodV0s();
1153
1154         gamma = new TClonesArray("AliKFParticle", gammaV0->GetEntriesFast());
1155         gamma->SetOwner(kTRUE);
1156
1157         // make a copy
1158         for(Int_t i=0; i < gammaV0->GetEntriesFast(); ++i)
1159         {
1160                 AliKFParticle* gamKF = (AliKFParticle*)gammaV0->At(i);
1161                 new ((*gamma)[i]) AliKFParticle(*gamKF);
1162                 posIndex.push_back(fV0Reader->GetPindex(i));
1163                 negIndex.push_back(fV0Reader->GetNindex(i));
1164         }
1165 }
1166
1167 //--------------------------------------------------------------------------
1168 TClonesArray* AliAnalysisTaskGammaConvDalitz::IndexToAliKFParticle(const vector<Int_t>& index, Int_t PDG)
1169 {
1170 //
1171 // Convert track index vector to AliKFParticle array
1172 //
1173         TClonesArray* indexKF = new TClonesArray("AliKFParticle",index.size());
1174         indexKF->SetOwner(kTRUE);
1175
1176         for(UInt_t i = 0; i < index.size(); ++i)
1177         {
1178                 AliESDtrack* t = fESDEvent->GetTrack(index[i]);
1179                 new((*indexKF)[i]) AliKFParticle(*t->GetConstrainedParam(), PDG);
1180         }
1181
1182         return indexKF;
1183 }
1184
1185 //--------------------------------------------------------------------------
1186 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindElectronFromPi0Dalitz(const vector<Int_t>& candidates, Int_t PDG)
1187 {
1188 //
1189 // Find true electrons from pi0 Dalitz decay candidates with MC
1190 //
1191         TClonesArray* elec = new TClonesArray("AliKFParticle");
1192         elec->SetOwner(kTRUE);
1193
1194         for(UInt_t i=0, j=0; i < candidates.size(); ++i)
1195         {
1196                 AliESDtrack* track = fESDEvent->GetTrack(candidates[i]);
1197                 Int_t trackLabel = TMath::Abs(track->GetLabel());
1198
1199                 if( fStack->Particle(trackLabel)->GetPdgCode() != PDG ) continue;
1200                 if( !IsPi0DalitzDaughter(trackLabel) ) continue;
1201
1202                 new ((*elec)[j++]) AliKFParticle(*track->GetConstrainedParam(), PDG);
1203         }
1204
1205         return elec;
1206 }
1207
1208 //--------------------------------------------------------------------------
1209 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindGammaFromPi0Dalitz(const TClonesArray* gamma, const vector<Int_t>& posIdx, const vector<Int_t>& negIdx)
1210 {
1211 //
1212 // Find true gammas from pi0 Dalitz decay candidates with MC
1213 //
1214         TClonesArray* gammaPi0 = new TClonesArray("AliKFParticle");
1215         gammaPi0->SetOwner(kTRUE);
1216
1217         for(Int_t i=0, j=0; i < gamma->GetEntriesFast(); ++i)
1218         {
1219                 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(i);
1220
1221                 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posIdx[i]))->GetLabel());
1222                 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negIdx[i]))->GetLabel());
1223
1224                 if( !HaveSameMother(labelv1,labelv2) ) continue;
1225
1226                 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
1227
1228                 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
1229
1230                 if( !IsPi0DalitzDaughter( labelGamma) ) continue;
1231
1232                 new ((*gammaPi0)[j++]) AliKFParticle(*gamKF);
1233         }
1234
1235         return gammaPi0;
1236 }
1237
1238 //--------------------------------------------------------------------------
1239 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindGamma(const TClonesArray* gamma, const vector<Int_t>& posIdx, const vector<Int_t>& negIdx)
1240 {
1241 //
1242 // Find true gammas from gamma candidates with MC
1243 //
1244         TClonesArray* gammaConv = new TClonesArray("AliKFParticle");
1245         gammaConv->SetOwner(kTRUE);
1246
1247         for(Int_t i=0, j=0; i < gamma->GetEntriesFast(); ++i)
1248         {
1249                 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(i);
1250
1251                 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posIdx[i]))->GetLabel());
1252                 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negIdx[i]))->GetLabel());
1253
1254                 if( !HaveSameMother(labelv1,labelv2) ) continue;
1255
1256                 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
1257
1258                 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
1259
1260                 new ((*gammaConv)[j++]) AliKFParticle(*gamKF);
1261         }
1262
1263         return gammaConv;
1264 }
1265
1266 //--------------------------------------------------------------
1267 void AliAnalysisTaskGammaConvDalitz::ESDtrackIndexCut(vector<Int_t>& pos, vector<Int_t>& neg, const TClonesArray* gamma)
1268 {
1269 //
1270 // Remove repeated electron candidate tracks
1271 // according to the gamma candidate array
1272 //
1273         vector<Bool_t> posTag(pos.size(),kTRUE);
1274         vector<Bool_t> negTag(neg.size(),kTRUE);
1275
1276         for(Int_t i=0; i < gamma->GetEntriesFast(); ++i)
1277         {
1278                 Int_t gamPosIndex = fGammaCandidatePosIndex[i];
1279                 Int_t gamNegIndex = fGammaCandidateNegIndex[i];
1280
1281                 for( UInt_t j=0; j < pos.size(); ++j )
1282                 {
1283                         if(pos[j] == gamPosIndex || pos[j] == gamNegIndex) posTag[j] = kFALSE;
1284                 }
1285
1286                 for( UInt_t j=0; j < neg.size(); ++j )
1287                 {
1288                         if(neg[j] == gamPosIndex || neg[j] == gamNegIndex) negTag[j] = kFALSE;
1289                 }
1290         }
1291
1292         CleanArray(pos, posTag);
1293         CleanArray(neg, negTag);
1294 }
1295
1296 //--------------------------------------------------------------------------
1297 void AliAnalysisTaskGammaConvDalitz::PsiPairCut(vector<Int_t>& pos, vector<Int_t>& neg)
1298 {
1299 //
1300 // Remove electron candidates from gamma conversions
1301 // according to the Psi pair angle
1302 //
1303     vector<Bool_t> posTag(pos.size(), kTRUE);
1304     vector<Bool_t> negTag(neg.size(), kTRUE);
1305
1306     for( UInt_t i=0; i < pos.size(); ++i )
1307     {
1308         AliESDtrack* posTrack = fESDEvent->GetTrack(pos[i]);
1309
1310         for( UInt_t j=0; j < neg.size(); ++j )
1311         {
1312             AliESDtrack* negTrack = fESDEvent->GetTrack(neg[j]);
1313
1314             Double_t psiPair = GetPsiPair(posTrack, negTrack);
1315             Double_t deltaPhi = fMagFieldSign * TVector2::Phi_mpi_pi( negTrack->GetConstrainedParam()->Phi()-posTrack->GetConstrainedParam()->Phi());
1316
1317             if(IsFromGammaConversion( psiPair, deltaPhi ))
1318             {
1319                 posTag[i] = kFALSE;
1320                 negTag[j] = kFALSE;
1321             }
1322         }
1323      }
1324
1325      CleanArray(pos, posTag);
1326      CleanArray(neg, negTag);
1327 }
1328
1329 //-----------------------------------------------------------------------------------
1330 void AliAnalysisTaskGammaConvDalitz::MassCut(vector<Int_t>& pos, vector<Int_t>& neg)
1331 {
1332 //
1333 // Remove electron candidates pairs
1334 // which have mass not in the range (fMassCutMin,fMassCutMax)
1335 //
1336     vector<Bool_t> posTag(pos.size(), kTRUE);
1337     vector<Bool_t> negTag(neg.size(), kTRUE);
1338
1339     for( UInt_t i=0; i < pos.size(); ++i )
1340     {
1341         AliESDtrack* posTrack = fESDEvent->GetTrack(pos[i]);
1342
1343         Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1344         TLorentzVector posLV;
1345         posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
1346
1347         for( UInt_t j=0; j < neg.size(); ++j )
1348         {
1349             AliESDtrack* negTrack = fESDEvent->GetTrack(neg[j]);
1350
1351             Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1352             TLorentzVector negLV;
1353             negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
1354
1355             TLorentzVector posnegLV = posLV + negLV;
1356
1357             if( (posnegLV.M() < fMassCutMin) || (posnegLV.M() > fMassCutMax) )
1358             {
1359                 posTag[i] = kFALSE;
1360                 negTag[j] = kFALSE;
1361             }
1362         }
1363      }
1364
1365      CleanArray(pos, posTag);
1366      CleanArray(neg, negTag);
1367 }
1368
1369 //-----------------------------------------------------------------------------------------------
1370 void AliAnalysisTaskGammaConvDalitz::CleanArray(vector<Int_t>& x, const vector<Bool_t>& tag)
1371 {
1372 //
1373 // Clean the x array according to the tag parameter
1374 //
1375         vector<Int_t> tmp;
1376
1377         for(UInt_t i=0; i< x.size(); ++i)
1378         {
1379                 if(tag[i]) tmp.push_back(x[i]);
1380         }
1381
1382         x = tmp;
1383 }
1384
1385 //--------------------------------------------------------------------------
1386 void AliAnalysisTaskGammaConvDalitz::AngleEposEnegGammaCut( const vector<Int_t>& posIdx, const vector<Int_t>& negIdx, const TClonesArray* candidates, TClonesArray*& gamma, vector<Int_t>& posGamIdx, vector<Int_t>& negGamIdx)
1387 {
1388 //
1389 // Remove gamma candidates according to
1390 // the angle between the plane e+,e- and the gamma
1391 //
1392         vector<Bool_t> gammaTag(candidates->GetEntriesFast(), kTRUE);
1393
1394         for( UInt_t iPos=0; iPos < posIdx.size(); ++iPos )
1395         {
1396                 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1397                 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1398                 TVector3 xMom(posMom[0],posMom[1],posMom[2]);
1399
1400                 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1401                 {
1402                         AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1403                         Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1404                         TVector3 yMom(negMom[0],negMom[1],negMom[2]);
1405
1406                         // normal vector to x+y- plane
1407                         TVector3 planePosNeg = xMom.Cross(yMom);
1408                         for(Int_t i=0; i < candidates->GetEntriesFast(); ++i)
1409                         {
1410                                 AliKFParticle* gamKF = (AliKFParticle*)candidates->At(i);
1411                                 TVector3 gamMom(gamKF->Px(),gamKF->Py(),gamKF->Pz());
1412                                 if (planePosNeg.Angle(gamMom) < 1. || planePosNeg.Angle(gamMom) > 2.)
1413                                 {
1414                                         gammaTag[i] = kFALSE;
1415                                 }
1416                         }
1417                 }
1418         }
1419
1420         // Rebuild gamma candidates array
1421
1422         if(gamma) delete gamma;
1423         gamma = new TClonesArray("AliKFParticle");
1424         gamma->SetOwner(kTRUE);
1425
1426         posGamIdx.clear();
1427         negGamIdx.clear();
1428
1429         for(Int_t i=0, j=0; i < candidates->GetEntriesFast(); ++i)
1430         {
1431                 AliKFParticle* iGamma = (AliKFParticle*)candidates->At(i);
1432                 if(gammaTag[i])
1433                 {
1434                         new ((*gamma)[j++]) AliKFParticle(*iGamma);
1435                         posGamIdx.push_back(fV0Reader->GetPindex(i));
1436                         negGamIdx.push_back(fV0Reader->GetNindex(i));
1437                 }
1438         }
1439 }
1440
1441 //--------------------------------------------------------------------------
1442 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindDalitzPair(const TClonesArray* pos, const TClonesArray* neg)
1443 {
1444 //
1445 // Find Dalitz pair candidates
1446 //
1447         TClonesArray* dalitz = new TClonesArray("TLorentzVector");
1448         dalitz->SetOwner(kTRUE);
1449
1450         for( Int_t iPos=0, j=0; iPos < pos->GetEntriesFast(); ++iPos )
1451         {
1452                 AliKFParticle* posKF = (AliKFParticle*)pos->At(iPos);
1453
1454                 TLorentzVector posLV;
1455                 posLV.SetXYZM(posKF->Px(),posKF->Py(),posKF->Pz(),fkElectronMass);
1456
1457                 for( Int_t iNeg=0; iNeg < neg->GetEntriesFast(); ++iNeg )
1458                 {
1459                         AliKFParticle* negKF = (AliKFParticle*)neg->At(iNeg);
1460
1461                         TLorentzVector negLV;
1462                         negLV.SetXYZM(negKF->Px(),negKF->Py(),negKF->Pz(),fkElectronMass);
1463
1464                         if(fUseAliKF)
1465                         {
1466                                 AliKFParticle posNeg( *posKF, *negKF);
1467
1468                                 TLorentzVector posNegLV;
1469                                 posNegLV.SetXYZM(posNeg.Px(), posNeg.Py(), posNeg.Pz(), posNeg.GetMass());
1470                                 new ((*dalitz)[j++]) TLorentzVector(posNegLV);
1471                         }
1472                         else
1473                         {
1474                                 new ((*dalitz)[j++]) TLorentzVector(posLV + negLV);
1475                         }
1476                 }
1477         }
1478
1479         return dalitz;
1480 }
1481
1482 //--------------------------------------------------------------------------
1483 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindPi0Dalitz(const TClonesArray* pos, const TClonesArray* neg, const TClonesArray* gamma)
1484 {
1485 //
1486 // Find pi0 Dalitz decay candidates
1487 //
1488         TClonesArray* pi0 = new TClonesArray("TLorentzVector");
1489         pi0->SetOwner(kTRUE);
1490
1491         for( Int_t iPos=0, j=0; iPos < pos->GetEntriesFast(); ++iPos )
1492         {
1493                 AliKFParticle* posKF = (AliKFParticle*)pos->At(iPos);
1494
1495                 TLorentzVector posLV;
1496                 posLV.SetXYZM(posKF->Px(),posKF->Py(),posKF->Pz(),fkElectronMass);
1497
1498                 for( Int_t iNeg=0; iNeg < neg->GetEntriesFast(); ++iNeg )
1499                 {
1500                         AliKFParticle* negKF = (AliKFParticle*)neg->At(iNeg);
1501
1502                         TLorentzVector negLV;
1503                         negLV.SetXYZM(negKF->Px(),negKF->Py(),negKF->Pz(),fkElectronMass);
1504
1505                         for(Int_t iGam=0; iGam < gamma->GetEntriesFast(); ++iGam)
1506                         {
1507                                 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(iGam);
1508
1509                                 if(fUseAliKF)
1510                                 {
1511                                         AliKFParticle posNegGam( *posKF, *negKF, *gamKF );
1512                                         TLorentzVector posNegGamLV;
1513                                         posNegGamLV.SetXYZM(posNegGam.Px(),posNegGam.Py(),posNegGam.Pz(),posNegGam.GetMass());
1514                                         new ((*pi0)[j++]) TLorentzVector(posNegGamLV);
1515                                 }
1516                                 else
1517                                 {
1518                                         TLorentzVector gamLV;
1519                                         gamLV.SetXYZM(gamKF->Px(),gamKF->Py(),gamKF->Pz(),0);
1520
1521                                         new ((*pi0)[j++]) TLorentzVector(posLV + negLV + gamLV);
1522                                 }
1523                         }
1524                 }
1525         }
1526
1527         return pi0;
1528 }
1529
1530 //--------------------------------------------------------------------------
1531 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindDalitzPair(const vector<Int_t>& posIdx, const vector<Int_t>& negIdx)
1532 {
1533 //
1534 // Find true Dalitz pairs from Dalitz pair candidats with MC
1535 //
1536         TClonesArray* dalitz = new TClonesArray("TLorentzVector");
1537         dalitz->SetOwner(kTRUE);
1538
1539         for( UInt_t iPos=0, j=0; iPos < posIdx.size(); ++iPos )
1540         {
1541                 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1542                 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1543                 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1544
1545                 TLorentzVector posLV;
1546                 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
1547
1548                 AliKFParticle posKF( *posTrack->GetConstrainedParam(), ::kPositron );
1549
1550                 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1551                 {
1552                         AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1553                         Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1554
1555                         if(!IsDalitzPair(posLabel,negLabel)) continue;
1556
1557                         if(fUseAliKF)
1558                         {
1559                                 AliKFParticle negKF( *negTrack->GetConstrainedParam(), ::kElectron );
1560                                 AliKFParticle posNeg( posKF, negKF);
1561
1562                                 TLorentzVector posNegLV;
1563                                 posNegLV.SetXYZM(posNeg.Px(),posNeg.Py(),posNeg.Pz(),posNeg.GetMass());
1564
1565                                 new ((*dalitz)[j++]) TLorentzVector(posNegLV);
1566                         }
1567                         else // TLorentzVector
1568                         {
1569                                 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1570
1571                                 TLorentzVector negLV;
1572                                 negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
1573
1574                                 new ((*dalitz)[j++]) TLorentzVector(posLV + negLV);
1575                         }
1576                 }
1577         }
1578
1579         return dalitz;
1580 }
1581
1582 //--------------------------------------------------------------------------
1583 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindPi0Dalitz(const vector<Int_t>& posIdx, const vector<Int_t>& negIdx, const TClonesArray* gamma, const vector<Int_t>& posGam, const vector<Int_t>& negGam)
1584 {
1585 //
1586 // Find true pi0 Dalitz decay from pi0 candidates with MC
1587 //
1588         TClonesArray* pi0 = new TClonesArray("TLorentzVector");
1589         pi0->SetOwner(kTRUE);
1590
1591         for( UInt_t iPos=0, j=0; iPos < posIdx.size(); ++iPos )
1592         {
1593                 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1594                 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1595                 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1596
1597                 TLorentzVector posLV;
1598                 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
1599
1600                 AliKFParticle posKF( *posTrack->GetConstrainedParam(), ::kPositron );
1601
1602                 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1603                 {
1604                         AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1605                         Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1606
1607                         if(!IsDalitzPair(posLabel,negLabel)) continue;
1608
1609                         Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1610
1611                         TLorentzVector negLV;
1612                         negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
1613
1614                         AliKFParticle negKF( *negTrack->GetConstrainedParam(), ::kElectron );
1615
1616                         for(Int_t iGam=0; iGam < gamma->GetEntriesFast(); ++iGam)
1617                         {
1618                                 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(iGam);
1619
1620                                 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posGam[iGam]))->GetLabel());
1621                                 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negGam[iGam]))->GetLabel());
1622
1623                                 if( !HaveSameMother(labelv1,labelv2) ) continue;
1624
1625                                 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
1626
1627                                 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
1628                                 if( !HaveSameMother(labelGamma, posLabel) ) continue;
1629
1630                                 if(fUseAliKF)
1631                                 {
1632                                         AliKFParticle posNegGam( posKF, negKF, *gamKF );
1633                                         TLorentzVector posNegGamLV;
1634                                         posNegGamLV.SetXYZM(posNegGam.Px(),posNegGam.Py(),posNegGam.Pz(),posNegGam.GetMass());
1635                                         new ((*pi0)[j++]) TLorentzVector(posNegGamLV);
1636                                 }
1637                                 else // TLorentzVector
1638                                 {
1639                                         TLorentzVector gamLV;
1640                                         gamLV.SetXYZM(gamKF->Px(),gamKF->Py(),gamKF->Pz(),0);
1641
1642                                         new ((*pi0)[j++]) TLorentzVector(posLV + negLV + gamLV);
1643                                 }
1644                         }
1645                 }
1646         }
1647
1648         return pi0;
1649 }
1650
1651 //-----------------------------------------------------------------------------------------------
1652 void AliAnalysisTaskGammaConvDalitz::UpdateGammaPool(const TClonesArray* gamma)
1653 {
1654 //
1655 // Update gamma event pool for background computation
1656 //
1657         if( fDebug ) AliInfo("=> UpdateGammaPool");
1658
1659         // cycle
1660         for(Int_t j=0; j< gamma->GetEntriesFast(); ++j)
1661         {
1662                 if((AliKFParticle*)fGammaPool->At(fGamPoolPos)) delete (AliKFParticle*)fGammaPool->RemoveAt(fGamPoolPos);
1663                 new ((*fGammaPool)[fGamPoolPos]) AliKFParticle( *((AliKFParticle*)gamma->At(j)));
1664                 ++fGamPoolPos;
1665                 if(fGamPoolPos == fPoolMaxSize)
1666                 {
1667                         fGamPoolPos = 0;
1668                 }
1669         }
1670 }
1671
1672 void AliAnalysisTaskGammaConvDalitz::UpdateElectronPool(TClonesArray* elec) // FIXME: const
1673 {
1674 //
1675 // Update electron event pool for background computation
1676 //
1677         Int_t multiplicity = fV0Reader->CountESDTracks();
1678
1679
1680         fBGEventHandler->AddElectronEvent(elec,fESDEvent->GetPrimaryVertex()->GetZ(),multiplicity);
1681 }
1682
1683 //-----------------------------------------------------------------------------------------------
1684 TClonesArray* AliAnalysisTaskGammaConvDalitz::GammasFromBGHandler() const
1685 {
1686 //
1687 // Gamma copy from events with same multiplicity and Z
1688 //
1689         if( fDebug ) AliInfo("=> GammasFromBGHandler");
1690
1691         Int_t zbin = fBGEventHandler->GetZBinIndex(fV0Reader->GetVertexZ());
1692         Int_t mbin = fBGEventHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
1693
1694
1695
1696         TClonesArray* gammaPool = new TClonesArray("AliKFParticle");
1697         gammaPool->SetOwner(kTRUE);
1698
1699         for( Int_t iEventBG=0; iEventBG < fV0Reader->GetNBGEvents(); ++iEventBG )
1700         {
1701                 AliGammaConversionKFVector* gammaV0s = fBGEventHandler->GetBGGoodV0s(zbin,mbin,iEventBG);
1702                 for( UInt_t i = 0; i < gammaV0s->size(); ++i)
1703                 {
1704                         new ((*gammaPool)[i]) AliKFParticle( *((AliKFParticle*)gammaV0s->at(i)) );
1705                 }
1706         }
1707
1708         return gammaPool;
1709 }
1710
1711 //-----------------------------------------------------------------------------------------------
1712 TClonesArray* AliAnalysisTaskGammaConvDalitz::ElectronFromBGHandler() const
1713 {
1714 //
1715 // Electron copy from events with same multiplicity and Z
1716 //
1717         if( fDebug ) AliInfo("=> ElectronFromBGHandler");
1718
1719         TClonesArray* electronPool = new TClonesArray("AliKFParticle");
1720         electronPool->SetOwner(kTRUE);
1721
1722         Int_t multiplicity = fV0Reader->CountESDTracks();
1723
1724
1725
1726
1727         for( Int_t iEventBG=0; iEventBG < fV0Reader->GetNBGEvents(); ++iEventBG )
1728         {
1729                 AliGammaConversionKFVector* electronNeg =  fBGEventHandler->GetBGGoodENeg(iEventBG,fESDEvent->GetPrimaryVertex()->GetZ(),multiplicity);
1730                 for (UInt_t i = 0; i < electronNeg->size(); ++i  )
1731                 {
1732                         new ((*electronPool)[i]) AliKFParticle( *((AliKFParticle*)electronNeg->at(i)) );
1733                 }
1734         }
1735
1736         return electronPool;
1737 }
1738
1739 //-----------------------------------------------------------------------------------------------
1740 Int_t AliAnalysisTaskGammaConvDalitz::GetMonteCarloPid(const AliESDtrack* t) const
1741 {
1742 //
1743 // Get track pid according to MC
1744 //
1745     Int_t label   = TMath::Abs(t->GetLabel());
1746     Int_t pdgCode = TMath::Abs(fStack->Particle(label)->GetPdgCode());
1747
1748     switch(pdgCode)
1749     {
1750          case ::kElectron:  return AliPID::kElectron;
1751          case ::kMuonMinus: return AliPID::kMuon;
1752          case ::kPiPlus:    return AliPID::kPion;
1753          case ::kKPlus:     return AliPID::kKaon;
1754          case ::kProton:    return AliPID::kProton;
1755     }
1756
1757     return -1;
1758 }
1759
1760 //-----------------------------------------------------------------------------------------------
1761 //FIXME PID ITS
1762 // NOTE prior should be estimated from data
1763 // NOTE: move to config
1764
1765 Int_t AliAnalysisTaskGammaConvDalitz::GetBayesPid(const AliESDtrack* t, Int_t trackType ) const
1766 {
1767 //
1768 // Get track pid according to Bayes' formula
1769 //
1770         double priors[AliPID::kSPECIES] = {0.009, 0.01, 0.82, 0.10, 0.05};
1771         Double_t detectoProb[AliPID::kSPECIES];
1772
1773         if( trackType == kITSsaTrack )  // ITS standalone pid
1774         {
1775            t->GetITSpid( detectoProb );
1776         }
1777         else  // global track
1778         {
1779            t->GetESDpid( detectoProb );
1780         }
1781
1782         AliPID bayesPID( detectoProb );
1783         return bayesPID.GetMostProbable( priors );
1784 }
1785
1786 //-----------------------------------------------------------------------------------------------
1787 Int_t AliAnalysisTaskGammaConvDalitz::GetNSigmaPid(const AliESDtrack* t, Int_t trackType ) const
1788 {
1789 //
1790 // Get track pid according to a n-sigma cut around ITS and/or TPC signals
1791 //
1792         if( trackType == kITSsaTrack)   // ITS standalone tracks
1793         {
1794             Double_t mom = t->GetP();
1795
1796
1797             // ITS Number of sigmas (FIXME: add new fESDpidCuts)
1798             // NOTE there is not AliESDpidCuts::SetITSnSigmaCut yet
1799             Double_t nElecSigma = fESDpid->NumberOfSigmasITS(t, AliPID::kElectron );
1800             Double_t nPionSigma = fESDpid->NumberOfSigmasITS(t, AliPID::kPion );
1801
1802             if( nElecSigma < 4. && nElecSigma > -3. && mom < .2  &&  nPionSigma < -1.5 )
1803             {
1804                 return AliPID::kElectron;
1805             }
1806         }
1807         else // global track
1808         {
1809             Double_t nElecSigma   = fESDpid->NumberOfSigmasTPC(t, AliPID::kElectron );
1810             Double_t nPionSigma   = fESDpid->NumberOfSigmasTPC(t, AliPID::kPion );
1811             Double_t nKaonSigma   = TMath::Abs(fESDpid->NumberOfSigmasTPC(t, AliPID::kKaon ));
1812             Double_t nProtonSigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(t, AliPID::kProton));
1813             if(     nElecSigma   > fNSigmaBelowElecTPCbethe && nElecSigma < fNSigmaAboveElecTPCbethe &&
1814                     nPionSigma   > fNSigmaAbovePionTPCbethe && //NOTE mom > 0.5
1815                     nKaonSigma   > fNSigmaAboveKaonTPCbethe &&
1816                     nProtonSigma > fNSigmaAboveProtonTPCbethe )
1817             {
1818                return AliPID::kElectron;
1819             }
1820             // NOTE: add other particle types
1821           }
1822
1823           return -1;
1824 }
1825
1826 //-----------------------------------------------------------------------------------------------
1827 Bool_t AliAnalysisTaskGammaConvDalitz::IsDalitzPair( Int_t posLabel, Int_t negLabel ) const
1828 {
1829 //
1830 // Returns true if the two particles is a Dalitz pair
1831 //
1832         if(!HaveSameMother(posLabel, negLabel)) return kFALSE;
1833
1834         TParticle* pos = fStack->Particle( posLabel );
1835         TParticle* neg = fStack->Particle( negLabel );
1836
1837         if( pos->GetPdgCode() != ::kPositron ) return kFALSE;
1838         if( neg->GetPdgCode() != ::kElectron ) return kFALSE;
1839
1840         //if( pos->GetUniqueID() != ::kPDecay ) return kFALSE;
1841         //if( neg->GetUniqueID() != ::kPDecay ) return kFALSE;
1842
1843         Int_t motherLabel = pos->GetMother(0);
1844         if( motherLabel < 0 ) return kFALSE;
1845
1846         TParticle* mother = fStack->Particle( motherLabel );
1847
1848         if( mother->GetPdgCode() != ::kPi0 ) return kFALSE;
1849         if( mother->GetNDaughters() != 3) return kFALSE;
1850         // NOTE: one of them must be a photon
1851
1852         return kTRUE;
1853 }
1854
1855 //-----------------------------------------------------------------------------------------------
1856 Bool_t AliAnalysisTaskGammaConvDalitz::IsPi0DalitzDaughter( Int_t label ) const
1857 {
1858 //
1859 // Returns true if the particle comes from Pi0 -> e+ e- gamma
1860 //
1861         Bool_t ePlusFlag  = kFALSE;
1862         Bool_t eMinusFlag = kFALSE;
1863         Bool_t gammaFlag  = kFALSE;
1864
1865         Int_t motherLabel = fStack->Particle( label )->GetMother(0);
1866
1867         if( motherLabel < 0 ) return kFALSE;
1868
1869         TParticle* mother = fStack->Particle( motherLabel );
1870
1871         if ( mother->GetPdgCode() != ::kPi0 ) return kFALSE;
1872
1873         if ( mother->GetNDaughters() != 3 ) return kFALSE;
1874
1875         for( Int_t idx = mother->GetFirstDaughter(); idx <= mother->GetLastDaughter(); ++idx )
1876         {
1877                 switch( fStack->Particle(idx)->GetPdgCode())
1878                 {
1879                         case ::kPositron:
1880                                 ePlusFlag  = kTRUE;
1881                                 break;
1882                         case ::kElectron:
1883                                 eMinusFlag = kTRUE;
1884                                 break;
1885                         case ::kGamma:
1886                                 gammaFlag  = kTRUE;
1887                                 break;
1888                 }
1889         }
1890
1891         return ( ePlusFlag && eMinusFlag && gammaFlag );
1892 }
1893
1894 //--------------------------------------------------------------------------
1895 Bool_t AliAnalysisTaskGammaConvDalitz::IsFromGammaConversion( Double_t psiPair, Double_t deltaPhi ) const
1896 {
1897 //
1898 // Returns true if it is a gamma conversion according to psi pair value
1899 //
1900         return ( (deltaPhi > fDeltaPhiCutMin  &&  deltaPhi < fDeltaPhiCutMax) &&
1901         TMath::Abs(psiPair) < ( fPsiPairCut - fPsiPairCut/fDeltaPhiCutMax * deltaPhi ) );
1902 }
1903
1904 //--------------------------------------------------------------------------
1905 Bool_t AliAnalysisTaskGammaConvDalitz::IsFromGammaConversion( Int_t posLabel, Int_t negLabel ) const
1906 {
1907 //
1908 // Returns true if it is a gamma conversion according to MC
1909 //
1910         if( !HaveSameMother(posLabel,negLabel) ) return kFALSE;
1911
1912         TParticle* pos = fStack->Particle( posLabel );
1913         TParticle* neg = fStack->Particle( negLabel );
1914
1915         if( pos->GetPdgCode() != ::kPositron ) return kFALSE;
1916         if( neg->GetPdgCode() != ::kElectron ) return kFALSE;
1917
1918         if( pos->GetUniqueID() != kPPair ) return kFALSE;
1919
1920         Int_t motherLabel = pos->GetMother(0);
1921         if( motherLabel < 0 ) return kFALSE;
1922
1923         TParticle* mother = fStack->Particle( motherLabel );
1924
1925         if( mother->GetPdgCode() != ::kGamma ) return kFALSE;
1926
1927         return kTRUE;
1928 }
1929
1930 //-----------------------------------------------------------------------------------------------
1931 Bool_t AliAnalysisTaskGammaConvDalitz::HaveSameMother( Int_t label1, Int_t label2 ) const
1932 {
1933 //
1934 // Returns true if the two particle have the same mother
1935 //
1936         if(fStack->Particle( label1 )->GetMother(0) < 0 ) return kFALSE;
1937         return (fStack->Particle( label1 )->GetMother(0) == fStack->Particle( label2 )->GetMother(0));
1938 }
1939
1940 //-----------------------------------------------------------------------------------------------
1941 Double_t AliAnalysisTaskGammaConvDalitz::GetPsiPair( const AliESDtrack* trackPos, const AliESDtrack* trackNeg ) const
1942 {
1943 //
1944 // This angle is a measure for the contribution of the opening in polar
1945 // direction Δ0 to the opening angle ξ Pair
1946 //
1947 // Ref. Measurement of photons via conversion pairs with the PHENIX experiment at RHIC
1948 //      Master Thesis. Thorsten Dahms. 2005
1949 // https://twiki.cern.ch/twiki/pub/ALICE/GammaPhysicsPublications/tdahms_thesis.pdf
1950 //
1951     Double_t momPos[3];
1952     Double_t momNeg[3];
1953     if( trackPos->GetConstrainedPxPyPz(momPos) == 0 ) trackPos->GetPxPyPz( momPos );
1954     if( trackNeg->GetConstrainedPxPyPz(momNeg) == 0 ) trackNeg->GetPxPyPz( momNeg );
1955
1956     TVector3 posDaughter;
1957     TVector3 negDaughter;
1958
1959     posDaughter.SetXYZ( momPos[0], momPos[1], momPos[2] );
1960     negDaughter.SetXYZ( momNeg[0], momNeg[1], momNeg[2] );
1961
1962     Double_t deltaTheta = negDaughter.Theta() - posDaughter.Theta();
1963     Double_t openingAngle =  posDaughter.Angle( negDaughter );  //TMath::ACos( posDaughter.Dot(negDaughter)/(negDaughter.Mag()*posDaughter.Mag()) );
1964     if( openingAngle < 1e-20 ) return 0.;
1965     Double_t psiAngle = TMath::ASin( deltaTheta/openingAngle );
1966
1967     return psiAngle;
1968 }
1969
1970 //-----------------------------------------------------------------------------------------------
1971 Double_t AliAnalysisTaskGammaConvDalitz::GetPsiPair(const AliKFParticle* xPos, const AliKFParticle* yNeg ) const
1972 {
1973 //
1974 // Get psi pair value
1975 //
1976     TVector3 pos(xPos->GetPx(), xPos->GetPy(), xPos->GetPz());
1977     TVector3 neg(yNeg->GetPx(), yNeg->GetPy(), yNeg->GetPz());
1978
1979     Double_t deltaTheta = neg.Theta() - pos.Theta();
1980     Double_t openingAngle = pos.Angle( neg );
1981
1982     if( openingAngle < 1e-20 ) return 0.;
1983
1984     return TMath::ASin( deltaTheta/openingAngle );
1985 }
1986
1987 //-----------------------------------------------------------------------------------------------
1988 Double_t AliAnalysisTaskGammaConvDalitz::GetPsiPair(const TLorentzVector* xPos, const TLorentzVector* yNeg ) const
1989 {
1990 //
1991 // Get psi pair value
1992 //
1993         Double_t deltaTheta = yNeg->Theta() - xPos->Theta();
1994         Double_t openingAngle = xPos->Angle( yNeg->Vect() );
1995
1996         if( openingAngle < 1e-20 ) return 0.;
1997
1998         return TMath::ASin( deltaTheta/openingAngle );;
1999 }
2000
2001 // tmp  NOTE: Should go to AliV0Reader
2002 //-----------------------------------------------------------------------------------------------
2003 /*
2004 Double_t AliAnalysisTaskGammaConvDalitz::GetPsiAngleV0(
2005             Double_t  radiusVO,                     // radius at XY of VO vertex
2006             const AliExternalTrackParam* trackPos,  // pos. track parm. at V0 vertex
2007             const AliExternalTrackParam* trackNeg  // neg. track parm. at V0 vertex
2008          )
2009 {
2010     // This angle is a measure for the contribution of the opening in polar
2011     // direction Δ0 to the opening angle ξ Pair
2012
2013     // Ref. Measurement of photons via conversion pairs with the PHENIX experiment at RHIC
2014     //      Master Thesis. Thorsten Dahms. 2005
2015     // https://twiki.cern.ch/twiki/pub/ALICE/GammaPhysicsPublications/tdahms_thesis.pdf
2016
2017     const Double_t kForceDeltaPhi = 1.2;
2018     static const Double_t kB2C = 0.299792458e-3; // Taken from AliVParticle
2019
2020     Double_t psiAngle = 0.;  // Default value
2021
2022     static Double_t MagnField = fESDEvent->GetMagneticField();
2023     static Double_t MagnFieldG = fESDEvent->GetMagneticField()*kB2C;
2024
2025     if( TMath::Abs( MagnField ) < 1e-20 ) return psiAngle; // Do nothing if 0 mag field
2026
2027     // compute propagation radius for a fixed angle
2028     Double_t Rpos = radiusVO + trackPos->Pt()*TMath::Sin( kForceDeltaPhi ) / MagnFieldG;
2029     Double_t Rneg = radiusVO + trackNeg->Pt()*TMath::Sin( kForceDeltaPhi ) / MagnFieldG;
2030
2031     Double_t MomPos[3];
2032     Double_t MomNeg[3];
2033     if( trackPos->GetPxPyPzAt( Rpos, MagnField, MomPos ) == 0 ) trackPos->GetPxPyPz( MomPos );
2034     if( trackNeg->GetPxPyPzAt( Rneg, MagnField, MomNeg ) == 0 ) trackNeg->GetPxPyPz( MomNeg );
2035
2036     TVector3 PosDaughter;
2037     TVector3 NegDaughter;
2038     PosDaughter.SetXYZ( MomPos[0], MomPos[1], MomPos[2] );
2039     NegDaughter.SetXYZ( MomNeg[0], MomNeg[1], MomNeg[2] );
2040
2041     Double_t deltaTheta = NegDaughter.Theta() - PosDaughter.Theta();
2042     Double_t chiPar = PosDaughter.Angle( NegDaughter );  //TMath::ACos( PosDaughter.Dot(NegDaughter) / (NegDaughter.Mag()*PosDaughter.Mag()) );
2043     if( chiPar < 1e-20 ) return psiAngle;
2044
2045     psiAngle = TMath::ASin( deltaTheta / chiPar );
2046
2047     return psiAngle;
2048 }
2049 */