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
 239
 240
 241     if( fDebug ) AliInfo("=> UserExec");
 242
 243     if( fV0Reader == 0 )
 244     {
 245         AliFatal("no pointer to AliV0Reader");
 246         return;
 247     }
 248
 249     // Create list of gamma candidates in standalone mode
 250     // otherwise use the created ones by AliAnalysisTaskGammaConversion
 251     if( fStandalone )
 252     {
 253
 254        AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 255        AliESDInputHandler *esdHandler=0;
 256         if ( (esdHandler=dynamic_cast<AliESDInputHandler*>(man->GetInputEventHandler())) && esdHandler->GetESDpid() ){
 257             AliV0Reader::SetESDpid(esdHandler->GetESDpid());
 258      } else {
 259     //load esd pid bethe bloch parameters depending on the existance of the MC handler
 260     // yes: MC parameters
 261     // no:  data parameters
 262         if (!AliV0Reader::GetESDpid()){
 263           if (MCEvent() ) {
 264                 AliV0Reader::InitESDpid();
 265             } else {
 266                 AliV0Reader::InitESDpid(1);
 267             }
 268         }
 269       }
 270
 271
 272         if (MCEvent() ) {
 273
 274     // To avoid crashes due to unzip errors. Sometimes the trees are not there.
 275         AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 276
 277         if (!mcHandler){
 278          AliError("Could not retrive MC event handler!");
 279          return;
 280         }
 281
 282         if (!mcHandler->InitOk() ){
 283
 284         return;
 285         }
 286         if (!mcHandler->TreeK() ){
 287         return;
 288         }
 289          if (!mcHandler->TreeTR() ) {
 290         return;
 291          }
 292       }
 293
 294
 295
 296
 297        fV0Reader->SetInputAndMCEvent( InputEvent(), MCEvent() );
 298        fV0Reader->Initialize();
 299     }
 300
 301     if( fV0Reader->CheckForPrimaryVertex() == kFALSE )
 302     {
 303         if( fDebug ) AliInfo("no contributors to primary vertex");
 304         return;
 305     }
 306
 307     if( fV0Reader->CheckForPrimaryVertexZ() == kFALSE  )
 308     {
 309
 310         if( fDebug ) AliInfo("z vertex out of range");
 311         return;
 312     }
 313
 314     // Get Pointers
 315    fBGEventHandler = fV0Reader->GetBGHandler();
 316    fESDpid = fV0Reader->GetESDpid();
 317    fESDEvent = fV0Reader->GetESDEvent();
 318    if(fDoMC && MCEvent())
 319    {
 320         fStack= MCEvent()->Stack();
 321         fGCMCEvent=MCEvent();
 322    }
 323
 324     // Read the magnetic field sign from ESD
 325     if ( fReadMagFieldSign == kTRUE )
 326     {
 327         fMagFieldSign = (fESDEvent->GetMagneticField() < 0) ? 1 : -1;
 328     }
 329
 330     // Process MC information
 331     if(fDoMC)
 332     {
 333         ProcessMCData();
 334     }
 335
 336     if( fStandalone ){
 337             while(fV0Reader->NextV0()){}; //SelectGammas
 338             fV0Reader->ResetV0IndexNumber();
 339     }
 340
 341     CreateListOfDalitzPairCandidates();
 342     ProcessGammaElectronsForDalitzAnalysis();
 343
 344     if ( fStandalone ){
 345
 346       fV0Reader->UpdateEventByEventData();
 347
 348     }
 349
 350     PostData( 1, fOutputContainer );
 351 }
 352
 353
 354 void AliAnalysisTaskGammaConvDalitz::Terminate(Option_t */*option*/)
 355 {
 356 //
 357     if( fDebug ) AliInfo("Not to do anything in Terminate");
 358 }
 359
 360 //-----------------------------------------------------------------------------------------------
 361 void AliAnalysisTaskGammaConvDalitz::AdoptITSsaTrackCuts( AliESDtrackCuts* esdCuts )
 362 {
 363 //
 364 // set user ITSsa track cuts
 365 //
 366     if( fITSsaTrackCuts ) delete fITSsaTrackCuts;
 367
 368     if( esdCuts )
 369     {
 370         fITSsaTrackCuts = esdCuts;
 371     }
 372     else
 373     {
 374         // default cuts
 375         fITSsaTrackCuts = new AliESDtrackCuts("Default ITSsa track cuts for Pi0 Dalitz decay");
 376
 377         fITSsaTrackCuts->SetEtaRange( -0.9, 0.9 );
 378         fITSsaTrackCuts->SetAcceptKinkDaughters(kFALSE);
 379
 380         fITSsaTrackCuts->SetMinNClustersITS(2);
 381         fITSsaTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
 382         fITSsaTrackCuts->SetRequireITSRefit(kTRUE);
 383
 384         fITSsaTrackCuts->SetRequireSigmaToVertex(kTRUE);
 385         fITSsaTrackCuts->SetMaxNsigmaToVertex(3);
 386
 387         fITSsaTrackCuts->SetRequireITSStandAlone(kTRUE);
 388     }
 389 }
 390
 391 //-----------------------------------------------------------------------------------------------
 392 void AliAnalysisTaskGammaConvDalitz::AdoptESDtrackCuts( AliESDtrackCuts* esdCuts )
 393 {
 394 //
 395 // set user global track cuts
 396 //
 397     if( fESDtrackCuts ) delete fESDtrackCuts;
 398
 399     if( esdCuts )
 400     {
 401         fESDtrackCuts = esdCuts;
 402     }
 403     else
 404     {
 405         //default cuts
 406         fESDtrackCuts = new AliESDtrackCuts("Default global track cuts for Pi0 Dalitz decay");
 407
 408         fESDtrackCuts->SetEtaRange( -0.9, 0.9 );
 409         fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
 410
 411         fESDtrackCuts->SetMinNClustersITS(2);
 412         fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kFirst);
 413         fESDtrackCuts->SetRequireITSRefit(kTRUE);
 414
 415         fESDtrackCuts->SetRequireSigmaToVertex(kTRUE);
 416         fESDtrackCuts->SetMaxNsigmaToVertex(3);
 417
 418         fESDtrackCuts->SetMinNClustersTPC(80);
 419         fESDtrackCuts->SetMaxChi2PerClusterTPC(4.);
 420         fESDtrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
 421         fESDtrackCuts->SetRequireTPCRefit(kTRUE);
 422     }
 423 }
 424
 425 //-----------------------------------------------------------------------------------------------
 426 void AliAnalysisTaskGammaConvDalitz::AdoptESDpidCuts( AliESDpidCuts* esdPIDCuts )
 427 {
 428 //
 429 // set user pid cuts
 430 //
 431     if( fESDpidCuts ) delete fESDpidCuts;
 432     if( esdPIDCuts )
 433     {
 434         fESDpidCuts = esdPIDCuts;
 435     }
 436     else // default cuts
 437     {
 438         fESDpidCuts = new AliESDpidCuts("Electrons", "Electron PID cuts");
 439      //   fESDpidCuts->SetTPCnSigmaCut(AliPID::kElectron, 3.);
 440         fESDpidCuts->SetTPCnSigmaCut(AliPID::kElectron, -4., 6.);
 441     }
 442 }
 443
 444 //-----------------------------------------------------------------------------------------------
 445 void AliAnalysisTaskGammaConvDalitz::ProcessMCData()
 446 {
 447 //
 448 // Process generation
 449 //
 450         if( fDebug ) AliInfo("=> ProcessMCData");
 451
 452         fHistograms->FillTable("Table_Generation", 0);  //number of events
 453
 454         for ( Int_t i = 0; i < fStack->GetNtrack(); i++ )
 455         {
 456                 TParticle* iParticle = fStack->Particle( i );
 457                 if( !iParticle ) continue;
 458
 459                 if ( i >= fStack->GetNprimary() ) continue;  // is primary?
 460                 if ( iParticle->GetPdgCode() != ::kPi0 ) continue;  // is Pi0?
 461
 462                 if( iParticle->GetNDaughters() == 2 &&
 463                     fStack->Particle(iParticle->GetFirstDaughter())->GetPdgCode() == ::kGamma &&
 464                     fStack->Particle(iParticle->GetLastDaughter())->GetPdgCode() == ::kGamma )
 465                 {
 466                         fHistograms->FillTable("Table_Generation", 1);  // pi0 -> gg
 467                 }
 468
 469                 if ( iParticle->GetNDaughters() != 3 ) continue;    // Num == 3 (e+,e-,gamma)
 470
 471                 // Check for Pi0 Dalitz decay
 472                 TParticle* eposPi0 = 0;
 473                 TParticle* enegPi0 = 0;
 474                 TParticle* gammaPi0 = 0;
 475
 476                 for( Int_t idxPi0 = iParticle->GetFirstDaughter(); idxPi0 <= iParticle->GetLastDaughter(); idxPi0++ )
 477                 {
 478                         switch(fStack->Particle(idxPi0)->GetPdgCode())
 479                         {
 480                                 case ::kPositron:
 481                                         eposPi0 = fStack->Particle(idxPi0);
 482                                         break;
 483                                 case ::kElectron:
 484                                         enegPi0 = fStack->Particle(idxPi0);
 485                                         break;
 486                                 case ::kGamma:
 487                                         gammaPi0 = fStack->Particle(idxPi0);
 488                                         break;
 489                         }
 490                 }
 491
 492                 if (eposPi0==0 || enegPi0==0 || gammaPi0==0) continue;
 493
 494                 // found a Pi0 Dalitz decay
 495
 496                 fHistograms->FillTable("Table_Generation", 2);
 497                 fHistograms->FillHistogram("MC_Pi0Dalitz_P", iParticle->P());
 498                 fHistograms->FillHistogram("MC_Pi0Dalitz_Pt", iParticle->Pt());
 499                 fHistograms->FillHistogram("MC_Pi0Dalitz_Eta", iParticle->Eta());
 500                 fHistograms->FillHistogram("MC_Pi0Dalitz_Pt_vs_Y", Rapidity(iParticle), iParticle->Pt());
 501                 fHistograms->FillHistogram("MC_EposDalitz_Pt", eposPi0->Pt());
 502                 fHistograms->FillHistogram("MC_EposDalitz_Eta", eposPi0->Eta());
 503                 fHistograms->FillHistogram("MC_EnegDalitz_Pt", enegPi0->Pt());
 504                 fHistograms->FillHistogram("MC_EnegDalitz_Eta", enegPi0->Eta());
 505                 fHistograms->FillHistogram("MC_GammaPi0Dalitz_Pt", gammaPi0->Pt());
 506                 fHistograms->FillHistogram("MC_GammaPi0Dalitz_Eta", gammaPi0->Eta());
 507
 508                 // Angle between the gamma and the plane e+e-
 509                 TVector3 ePosMom( eposPi0->Px(), eposPi0->Py(), eposPi0->Pz() );
 510                 TVector3 eNegMom( enegPi0->Px(), enegPi0->Py(), enegPi0->Pz() );
 511                 TVector3 gamMom( gammaPi0->Px(), gammaPi0->Py() , gammaPi0->Pz() );
 512                 TVector3 planeEposEneg =  eNegMom.Cross( ePosMom );
 513                 Double_t anglePlaneGamma = planeEposEneg.Angle(gamMom);
 514
 515                 fHistograms->FillHistogram("MC_EposEnegDalitz_Angle", ePosMom.Angle(eNegMom) );
 516
 517                 fHistograms->FillHistogram("MC_EposEnegDalitz_GammaPi0_Angle", anglePlaneGamma);
 518                 fHistograms->FillHistogram("MC_EposEnegDalitz_GammaPi0_Angle_vs_P", anglePlaneGamma, gammaPi0->P());
 519                 fHistograms->FillHistogram("MC_EposEnegDalitz_GammaPi0_Angle_vs_Pt", anglePlaneGamma, gammaPi0->Pt());
 520
 521
 522         // check for gamma conversion
 523         Bool_t daugGammaElectron    = kFALSE;
 524         Bool_t daugGammaPositron    = kFALSE;  // acceptance
 525         Bool_t daugGammaElectronAll = kFALSE;
 526         Bool_t daugGammaPositronAll = kFALSE;
 527
 528         // is the gamma converted? -> has 2 daughter e+e-
 529         // are e+ e- from gamma in the acceptance for the V0s
 530         if( gammaPi0->GetNDaughters() >= 2 )
 531         {
 532             for( Int_t tIndex=gammaPi0->GetFirstDaughter(); tIndex<=gammaPi0->GetLastDaughter(); ++tIndex )
 533             {
 534                 TParticle* tmpDaughter = fStack->Particle(tIndex);
 535
 536                 if( tmpDaughter->GetUniqueID() != kPPair ) continue; // check if the daughters come from a conversion
 537                 if( tmpDaughter->GetPdgCode() == ::kElectron )
 538                 { // e+
 539                     daugGammaElectronAll = kTRUE;
 540
 541                     if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() &&
 542                         ((TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()) < tmpDaughter->R() &&
 543                         (tmpDaughter->R()< fV0Reader->GetMaxRCut() ) )
 544                     {
 545                         daugGammaElectron = kTRUE;
 546                     }
 547                 }
 548                 else if( tmpDaughter->GetPdgCode() == ::kPositron )
 549                 {
 550                      daugGammaPositronAll = kTRUE;
 551                     if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() &&
 552                         ((TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()) < tmpDaughter->R() &&
 553                         (tmpDaughter->R() < fV0Reader->GetMaxRCut() ) )
 554                     {
 555                         daugGammaPositron = kTRUE;
 556                     }
 557                 }
 558             }
 559         }
 560
 561
 562        if(  daugGammaElectronAll && daugGammaPositronAll )
 563        {
 564           TParticle* tmpDaughter = fStack->Particle( gammaPi0->GetFirstDaughter() );
 565           fHistograms->FillHistogram("MC_GC_GammaPi0Dalitz_All_Z_vs_R",tmpDaughter->Vz(),tmpDaughter->R() );
 566        }
 567
 568         Float_t  etaMin, etaMax;
 569         fESDtrackCuts->GetEtaRange( etaMin, etaMax );
 570
 571         // e+e- pair within acceptance
 572         if ( TMath::Abs( eposPi0->Eta() ) < etaMax  && TMath::Abs( enegPi0->Eta() ) < etaMax )
 573         {
 574             fHistograms->FillHistogram("MC_Acceptance_EposDalitz_Pt", eposPi0->Pt());
 575             fHistograms->FillHistogram("MC_Acceptance_EposDalitz_Eta", eposPi0->Eta());
 576             fHistograms->FillHistogram("MC_Acceptance_EnegDalitz_Pt", enegPi0->Pt());
 577             fHistograms->FillHistogram("MC_Acceptance_EnegDalitz_Eta", enegPi0->Eta());
 578             fHistograms->FillHistogram("MC_Acceptance_DalitzPair_EposPt_vs_EnegPt", eposPi0->Pt(), enegPi0->Pt());
 579         }
 580
 581         // Pi0 (e+e-gamma) within acceptance
 582         //cout<<"Gamma Eta Cut"<<fV0Reader->GetEtaCut()<<endl;
 583
 584         if ( ( TMath::Abs( gammaPi0->Eta() ) < fV0Reader->GetEtaCut() && gammaPi0->R() < fV0Reader->GetMaxRCut() ) &&
 585              TMath::Abs( eposPi0->Eta() ) < etaMax  && TMath::Abs( enegPi0->Eta() ) < etaMax )
 586         {
 587             fHistograms->FillTable("Table_Generation",3);  //
 588             fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Pt",iParticle->Pt());
 589             fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Eta",iParticle->Eta());
 590             fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Radius",iParticle->R());
 591             fHistograms->FillHistogram("MC_Acceptance_GammaPi0Dalitz_Pt",gammaPi0->Pt());
 592             fHistograms->FillHistogram("MC_Acceptance_GammaPi0Dalitz_Eta",gammaPi0->Eta());
 593             fHistograms->FillHistogram("MC_Acceptance_EposPi0Dalitz_Pt",eposPi0->Pt());
 594             fHistograms->FillHistogram("MC_Acceptance_EposPi0Dalitz_Eta",eposPi0->Eta());
 595             fHistograms->FillHistogram("MC_Acceptance_EnegPi0Dalitz_Pt",enegPi0->Pt());
 596             fHistograms->FillHistogram("MC_Acceptance_EnegPi0Dalitz_Eta",enegPi0->Eta());
 597             fHistograms->FillHistogram("MC_Acceptance_DalitzPair_OpeningAngle", ePosMom.Angle(eNegMom) );
 598             fHistograms->FillHistogram("MC_Acceptance_Pi0Dalitz_Pt_vs_Y", Rapidity(iParticle), iParticle->Pt());
 599
 600            // Pi0 within acceptance with gamma converted
 601
 602             if ( daugGammaElectron && daugGammaPositron )
 603             {
 604                 fHistograms->FillTable("Table_Generation",4); //
 605
 606                 fHistograms->FillHistogram("MC_Acceptance_GC_Pi0Dalitz_Pt",iParticle->Pt());
 607                 fHistograms->FillHistogram("MC_Acceptance_GC_Pi0Dalitz_Eta",iParticle->Eta());
 608                 fHistograms->FillHistogram("MC_Acceptance_GC_EposPi0Dalitz_Pt",eposPi0->Pt());
 609                 fHistograms->FillHistogram("MC_Acceptance_GC_EposPi0Dalitz_Eta",eposPi0->Eta());
 610                 fHistograms->FillHistogram("MC_Acceptance_GC_EnegPi0Dalitz_Pt",enegPi0->Pt());
 611                 fHistograms->FillHistogram("MC_Acceptance_GC_EnegPi0Dalitz_Eta",enegPi0->Eta());
 612                 fHistograms->FillHistogram("MC_Acceptance_GC_GammaPi0Dalitz_Pt",gammaPi0->Pt());
 613                 fHistograms->FillHistogram("MC_Acceptance_GC_GammaPi0Dalitz_Eta",gammaPi0->Eta());
 614                 //fHistograms->FillHistogram("MC_Acceptance_GC_Gamma_Angle",anglePlaneGamma);
 615                 //fHistograms->FillHistogram("MC_Acceptance_GC_Gamma_Angle_vs_Pt",anglePlaneGamma,gammaPi0->Pt());
 616                 TParticle* tmpDaughter = fStack->Particle( gammaPi0->GetFirstDaughter() );
 617                 fHistograms->FillHistogram("MC_Acceptance_GC_GammaPi0Dalitz_Z_vs_R",tmpDaughter->Vz(),tmpDaughter->R() );
 618                 fHistograms->FillHistogram("MC_Acceptance_GC_Pi0Dalitz_Pt_vs_Y", Rapidity(iParticle), iParticle->Pt());
 619             }
 620         }
 621     }
 622 }
 623
 624 //-----------------------------------------------------------------------------------------------
 625 void AliAnalysisTaskGammaConvDalitz::CreateListOfDalitzPairCandidates()
 626 {
 627 //
 628 // Dalitz pair candidates
 629 //
 630         if( fDebug ) AliInfo("=> CreateListOfDalitzPairCandidates");
 631
 632         fEposCandidateIndex.clear();
 633         fEnegCandidateIndex.clear();
 634
 635         fHistograms->FillTable("Table_Cuts", 0);
 636
 637         for( Int_t i = 0; i < fESDEvent->GetNumberOfTracks(); ++i )
 638         {
 639                 AliESDtrack* iTrack = fESDEvent->GetTrack(i);
 640                 if ( !iTrack ) continue;
 641
 642
 643                 Double_t p[3];
 644
 645                 if ( !iTrack->GetConstrainedPxPyPz(p) ) continue;
 646
 647                 TVector3 iMom(p[0],p[1],p[2]);
 648
 649                 //
 650                 // Check track cuts and find track type
 651                 //
 652
 653                 Bool_t isTrackAccepted = 0;
 654                 Int_t trackType = -1;
 655                 switch(fTrkSelectionCriteria)
 656                 {
 657                         case kITSsaTrack:
 658                                 isTrackAccepted = fITSsaTrackCuts->AcceptTrack( iTrack );
 659                                 trackType = kITSsaTrack;
 660                                 break;
 661
 662                         case kGlobalTrack:
 663                                 isTrackAccepted = fESDtrackCuts->AcceptTrack( iTrack );
 664                                 trackType = kGlobalTrack;
 665                                 break;
 666
 667                         case kITSsaGlobalTrack:
 668                                 if(fITSsaTrackCuts->AcceptTrack( iTrack ) || fESDtrackCuts->AcceptTrack( iTrack ))
 669                                 {
 670                                         isTrackAccepted = kTRUE;
 671                                         if(fITSsaTrackCuts->AcceptTrack( iTrack )) trackType = kITSsaTrack;
 672                                         else trackType = kGlobalTrack;
 673                                 }
 674                                 break;
 675                 }
 676
 677                 if(!isTrackAccepted) continue;
 678
 679                 //
 680                 // PID
 681                 //
 682
 683                 Int_t pid=-1;
 684                 Int_t pidMC=-1;
 685
 686                 if(fUseBayesPID)
 687                 {
 688                         pid = GetBayesPid(iTrack,trackType);
 689                 }
 690                 else
 691                 {
 692                         pid = GetNSigmaPid(iTrack,trackType);
 693                 }
 694
 695                 if( fDoMC )
 696                 {
 697                         pidMC = GetMonteCarloPid(iTrack);
 698                         // pid table
 699                         Int_t iLabel = TMath::Abs(iTrack->GetLabel());
 700                         TParticle* iParticle = fStack->Particle(iLabel);
 701                         FillPidTable(iParticle, pid);
 702                 }
 703
 704                 // ITS standalone tracks
 705                 if( trackType == kITSsaTrack)
 706                 {
 707                         Double_t mom = iTrack->GetP();
 708                         Double_t signal = iTrack->GetITSsignal();
 709
 710                         fHistograms->FillHistogram( "ESD_ITSsa_dEdx_vs_P", mom, signal );
 711
 712                         if( pid == AliPID::kElectron )
 713                         {
 714
 715                                 fHistograms->FillHistogram( "ESD_ITSsa_PidCut_dEdx_vs_P", mom, signal );
 716                                 if(fDoMC && pid == pidMC)
 717                                 {
 718                                         fHistograms->FillHistogram( "MC_ESD_ITSsa_PidCut_dEdx_vs_P", mom, signal );
 719                                 }
 720                         }
 721
 722                         if( fDoMC && pidMC == AliPID::kElectron)
 723                         {
 724                                 fHistograms->FillHistogram( "MC_ESD_ITSsa_Electron_dEdx_vs_P", mom, signal );
 725                         }
 726                 }
 727
 728                 else  // global tracks
 729                 {
 730                         const AliExternalTrackParam *in = iTrack->GetInnerParam();
 731                         Double_t mom = in->GetP();
 732                         Double_t signal = iTrack->GetTPCsignal();
 733
 734                         fHistograms->FillHistogram( "ESD_TPC_dEdx_vs_P", mom, signal );
 735
 736                         if( fDoMC && pidMC == AliPID::kElectron )
 737                         {
 738                                 fHistograms->FillHistogram( "MC_ESD_TPC_Electron_dEdx_vs_P", mom, signal );
 739                         }
 740
 741                         if( pid == AliPID::kElectron )
 742                         {
 743                                 fHistograms->FillHistogram( "ESD_TPC_PidCut_dEdx_vs_P", mom, signal );
 744                                 if(fDoMC && pid == pidMC)
 745                                 {
 746                                         fHistograms->FillHistogram( "MC_ESD_TPC_PidCut_dEdx_vs_P", mom, signal );
 747                                 }
 748                         }
 749                 }
 750
 751                 if( AliPID::kElectron != pid) continue;
 752
 753                 // electron track candidates from here
 754
 755                 if( iTrack->GetSign() > 0 )
 756                 {
 757                         fEposCandidateIndex.push_back(i);
 758                 }
 759                 else
 760                 {
 761                         fEnegCandidateIndex.push_back(i);
 762                 }
 763         }
 764
 765         // gamma candidates
 766         GetGammaCandidates(fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
 767
 768         if(fDoMC)
 769         {
 770                 TClonesArray* pi0Dalitz = FindParticleDalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,1);
 771                 ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(1.,(Double_t)pi0Dalitz->GetEntriesFast());
 772                 delete pi0Dalitz;
 773         }
 774
 775         if(fUseTrackIndexCut) // remove repeated tracks
 776         {
 777                 ESDtrackIndexCut(fEposCandidateIndex,fEnegCandidateIndex, fGammaCandidates);
 778
 779                 if(fDoMC)
 780                 {
 781                         TClonesArray* pi0Dalitz = FindParticleDalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,1);
 782                         ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(2.,(Double_t)pi0Dalitz->GetEntriesFast());
 783                         delete pi0Dalitz;
 784                 }
 785         }
 786
 787         if(fUsePsiPairCut) // remove electrons from gamma conversions
 788         {
 789                 PsiPairCut(fEposCandidateIndex,fEnegCandidateIndex);
 790
 791                 if(fDoMC)
 792                 {
 793                         TClonesArray* pi0Dalitz = FindParticleDalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,1);
 794                         ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(3.,(Double_t)pi0Dalitz->GetEntriesFast());
 795                         delete pi0Dalitz;
 796                 }
 797         }
 798
 799         if( fUseMassCut )
 800         {
 801                 MassCut(fEposCandidateIndex, fEnegCandidateIndex);
 802
 803                 if(fDoMC)
 804                 {
 805                         TClonesArray* pi0Dalitz = FindParticleDalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,1);
 806                         ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(4.,(Double_t)pi0Dalitz->GetEntriesFast());
 807                         delete pi0Dalitz;
 808                 }
 809         }
 810
 811         if(fUseGammaCut)
 812         {
 813                 AngleEposEnegGammaCut(fEposCandidateIndex,fEnegCandidateIndex,fV0Reader->GetCurrentEventGoodV0s(), fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
 814
 815                 if(fDoMC)
 816                 {
 817                         TClonesArray* pi0Dalitz = FindParticleDalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,1);
 818                         ((TH1*)fHistograms->GetValue("Table_Cuts"))->Fill(5.,(Double_t)pi0Dalitz->GetEntriesFast());
 819                         delete pi0Dalitz;
 820                 }
 821         }
 822 }
 823
 824 //-----------------------------------------------------------------------------------------------
 825 void AliAnalysisTaskGammaConvDalitz::ProcessGammaElectronsForDalitzAnalysis()
 826 {
 827 //
 828 // Process gamma and electrons for pi0 Dalitz decay
 829 //
 830         if( fDebug ) AliInfo("=> ProcessGammaElectronsForDalitzAnalysis");
 831
 832
 833         fHistograms->FillTable( "Table_Reconstruction", 0); // number of events
 834
 835
 836
 837         TClonesArray* ePosCandidates = IndexToAliKFParticle(fEposCandidateIndex, ::kPositron);
 838
 839         for(Int_t i=0; i < ePosCandidates->GetEntriesFast(); ++i)
 840         {
 841                 AliKFParticle* epos = (AliKFParticle*) ePosCandidates->At(i);
 842                 fHistograms->FillHistogram("ESD_EposCandidates_Pt", epos->GetPt());
 843                 fHistograms->FillHistogram("ESD_EposCandidates_Eta", epos->GetEta());
 844                 fHistograms->FillTable( "Table_Reconstruction", 1);
 845         }
 846
 847         TClonesArray* eNegCandidates = IndexToAliKFParticle(fEnegCandidateIndex, ::kElectron);
 848
 849         for(Int_t i=0; i < eNegCandidates->GetEntriesFast(); ++i)
 850         {
 851                 AliKFParticle* eneg = (AliKFParticle*) eNegCandidates->At(i);
 852                 fHistograms->FillHistogram("ESD_EnegCandidates_Pt", eneg->GetPt());
 853                 fHistograms->FillHistogram("ESD_EnegCandidates_Eta", eneg->GetEta());
 854                 fHistograms->FillTable( "Table_Reconstruction", 2);
 855         }
 856
 857         TClonesArray* dalitzPairCandidates = FindDalitzPair(ePosCandidates, eNegCandidates);
 858         for(Int_t i=0; i < dalitzPairCandidates->GetEntriesFast(); ++i)
 859         {
 860                 TLorentzVector* dalitz = (TLorentzVector*)dalitzPairCandidates->At(i);
 861
 862                 fHistograms->FillHistogram("ESD_DalitzPairCandidates_Pt", dalitz->Pt());
 863                 fHistograms->FillHistogram("ESD_DalitzPairCandidates_InvMass", dalitz->M());
 864                 fHistograms->FillHistogram("ESD_DalitzPairCandidates_InvMass_vs_Pt",dalitz->M(),dalitz->Pt());
 865         }
 866
 867         // gamma candidates
 868         for(Int_t i=0; i < fGammaCandidates->GetEntriesFast(); ++i)
 869         {
 870                 AliKFParticle* gamma = (AliKFParticle*) fGammaCandidates->At(i);
 871                 fHistograms->FillHistogram("ESD_GammaCandidates_Pt", gamma->GetPt());
 872                 fHistograms->FillHistogram("ESD_GammaCandidates_Eta", gamma->GetEta());
 873         }
 874
 875         // psi pair for all candidates
 876         //if(fUsePsiPairCut)
 877         FillPsiPair(ePosCandidates,eNegCandidates,"ESD_EposEneg_PsiPair_vs_DPhi");
 878
 879         // Angle epos,eneg gamma
 880         FillAngle(ePosCandidates, fGammaCandidates, "ESD_EposEneg_GammaCandidates_Angle");
 881         FillAngle(eNegCandidates, fGammaCandidates, "ESD_EposEneg_GammaCandidates_Angle");
 882
 883         TClonesArray* pi0Candidates = FindParticleDalitz(ePosCandidates, eNegCandidates, fGammaCandidates,0);
 884         for(Int_t i=0; i < pi0Candidates->GetEntriesFast(); ++i)
 885         {
 886                 TLorentzVector* pi0 = (TLorentzVector*)pi0Candidates->At(i);
 887
 888                 fHistograms->FillHistogram("ESD_Pi0_P", pi0->P());
 889                 fHistograms->FillHistogram("ESD_Pi0_Pt", pi0->Pt());
 890                 fHistograms->FillHistogram("ESD_Pi0_Eta", pi0->Eta());
 891                 fHistograms->FillHistogram("ESD_Pi0_Y", pi0->Rapidity());
 892                 fHistograms->FillHistogram("ESD_Pi0_Phi", pi0->Phi());
 893                 fHistograms->FillHistogram("ESD_Pi0_Pt_vs_Y",pi0->Pt(),pi0->Rapidity());
 894                 fHistograms->FillHistogram("ESD_Pi0_InvMass", pi0->M());
 895                 fHistograms->FillHistogram("ESD_Pi0_InvMass_vs_Pt", pi0->M(),pi0->Pt());
 896                 fHistograms->FillHistogram("ESD_Pi0_InvMass_vs_Y",pi0->M(),pi0->Rapidity());
 897                 fHistograms->FillHistogram("ESD_Pi0_InvMass_vs_Eta",pi0->M(),pi0->Eta());
 898         }
 899
 900         for(Int_t iPos=0; iPos < ePosCandidates->GetEntriesFast(); ++iPos)
 901         {
 902                 AliKFParticle* lPosKF = (AliKFParticle*)ePosCandidates->At(iPos);
 903
 904                 for(Int_t iNeg=0; iNeg < eNegCandidates->GetEntriesFast(); ++iNeg)
 905                 {
 906                     AliKFParticle* lNegKF = (AliKFParticle*)eNegCandidates->At(iNeg);
 907                     AliKFParticle lPosNeg(*lPosKF,*lNegKF );
 908
 909                     for(Int_t iGam=0; iGam < fGammaCandidates->GetEntriesFast(); ++iGam)
 910                     {
 911                         AliKFParticle* lGamKF = (AliKFParticle*)fGammaCandidates->At(iGam);
 912
 913                         AliKFParticle lPosNegGam( *lPosKF, *lNegKF, *lGamKF );
 914
 915                         Double_t lDiffMass = lPosNegGam.GetMass() - lPosNeg.GetMass();
 916
 917                         fHistograms->FillHistogram("ESD_EposEnegGamma_InvMass_Diff",lDiffMass );
 918                         fHistograms->FillHistogram("ESD_EposEnegGamma_InvMass_vs_Pt_Diff",lDiffMass,lPosNegGam.GetPt());
 919                         fHistograms->FillHistogram("ESD_EposEnegGamma_InvMass",lPosNegGam.GetMass());
 920                         fHistograms->FillHistogram("ESD_EposEnegGamma_InvMass_vs_Pt",lPosNegGam.GetMass(),lPosNegGam.GetPt());
 921
 922                     }
 923                }
 924         }
 925
 926
 927
 928
 929         delete dalitzPairCandidates;
 930         delete pi0Candidates;
 931
 932         if(fComputeBkg)
 933         {
 934
 935                 // 1) e+e- dalitz
 936                 for(Int_t i=0; i < ePosCandidates->GetEntriesFast(); ++i)
 937                 {
 938                         AliKFParticle* epos1 = (AliKFParticle*) ePosCandidates->At(i);
 939
 940                     for(Int_t j=i+1; j < ePosCandidates->GetEntriesFast(); ++j)
 941                     {
 942                         AliKFParticle* epos2 = (AliKFParticle*) ePosCandidates->At(j);
 943                         AliKFParticle ePosePos( *epos1,*epos2 );
 944                         fHistograms->FillHistogram("ESD_BKG_LikeSign_InvMass",ePosePos.GetMass());
 945                         fHistograms->FillHistogram("ESD_BKG_LikeSign_InvMass_vs_Pt",ePosePos.GetMass(),ePosePos.GetPt());
 946
 947
 948                     }
 949
 950
 951                 }
 952                 for(Int_t i=0; i < eNegCandidates->GetEntriesFast(); ++i)
 953                 {
 954                         AliKFParticle* eneg1 = (AliKFParticle*) eNegCandidates->At(i);
 955
 956                     for(Int_t j=i+1; j < eNegCandidates->GetEntriesFast(); ++j)
 957                     {
 958                         AliKFParticle* eneg2 = (AliKFParticle*) eNegCandidates->At(j);
 959                         AliKFParticle eNegeNeg( *eneg1,*eneg2 );
 960                         fHistograms->FillHistogram("ESD_BKG_LikeSign_InvMass",eNegeNeg.GetMass());
 961                         fHistograms->FillHistogram("ESD_BKG_LikeSign_InvMass_vs_Pt",eNegeNeg.GetMass(),eNegeNeg.GetPt());
 962                     }
 963
 964
 965                 }
 966
 967                 // 1) e+e- with with gammas used in the signal
 968                 TClonesArray* pi0Bkg = FindParticleDalitz(ePosCandidates, eNegCandidates, fGammaPool,1);
 969
 970                 for(Int_t i=0; i < pi0Bkg->GetEntriesFast(); ++i)
 971                 {
 972                         TLorentzVector* pi0 = (TLorentzVector*)pi0Bkg->At(i);
 973                         fHistograms->FillHistogram("ESD_BKG_PrevGamma_InvMass", pi0->M());
 974                         fHistograms->FillHistogram("ESD_BKG_PrevGamma_InvMass_vs_Pt",pi0->M(),pi0->Pt());
 975                 }
 976                 ///////////////////////////////Temporal for Dalitz
 977
 978
 979
 980
 981
 982
 983
 984                 if(ePosCandidates->GetEntriesFast() > 0 &&
 985                    eNegCandidates->GetEntriesFast() > 0 &&
 986                    fGammaCandidates->GetEntriesFast() > 0)
 987                 {
 988                         UpdateGammaPool(fGammaCandidates);
 989                 }
 990
 991                 delete pi0Bkg;
 992
 993                 // 2) e+e- with gammas from a pool of events
 994                 TClonesArray* gammaBGHandler = GammasFromBGHandler();
 995                 pi0Bkg = FindParticleDalitz(ePosCandidates, eNegCandidates, gammaBGHandler,2);
 996
 997                 for(Int_t i=0; i < pi0Bkg->GetEntriesFast(); ++i)
 998                 {
 999                         TLorentzVector* pi0 = (TLorentzVector*)pi0Bkg->At(i);
1000                         fHistograms->FillHistogram("ESD_BKG_BGHandler_InvMass", pi0->M());
1001                         fHistograms->FillHistogram("ESD_BKG_BGHandler_InvMass_vs_Pt",pi0->M(), pi0->Pt());
1002                 }
1003
1004                 delete pi0Bkg;
1005
1006                 // 3) e+ with e-, gamma from a pool of events
1007                 TClonesArray* elecBGHandler = ElectronFromBGHandler();
1008                 pi0Bkg = FindParticleDalitz(ePosCandidates, elecBGHandler, gammaBGHandler,3);
1009
1010                 for(Int_t i=0; i < pi0Bkg->GetEntriesFast(); ++i)
1011                 {
1012                         TLorentzVector* pi0 = (TLorentzVector*)pi0Bkg->At(i);
1013                         fHistograms->FillHistogram("ESD_BKG_Electron_InvMass", pi0->M());
1014                         fHistograms->FillHistogram("ESD_BKG_Electron_InvMass_vs_Pt",pi0->M(), pi0->Pt());
1015                 }
1016
1017                 if(eNegCandidates->GetEntriesFast() > 0)
1018                 {
1019                         UpdateElectronPool(eNegCandidates);
1020                 }
1021
1022                 delete gammaBGHandler;
1023                 delete elecBGHandler;
1024                 delete pi0Bkg;
1025
1026         }
1027
1028         delete ePosCandidates;
1029         delete eNegCandidates;
1030
1031         if(fDoMC)
1032         {
1033                 TClonesArray* ePosPi0Dalitz = FindElectronFromPi0Dalitz(fEposCandidateIndex, ::kPositron);
1034                 for(Int_t i=0; i < ePosPi0Dalitz->GetEntriesFast(); ++i)
1035                 {
1036                         AliKFParticle* epos = (AliKFParticle*) ePosPi0Dalitz->At(i);
1037                         fHistograms->FillHistogram("MC_ESD_Pi0_EposDalitz_Pt", epos->GetPt());
1038                         fHistograms->FillHistogram("MC_ESD_Pi0_EposDalitz_Eta", epos->GetEta());
1039                         fHistograms->FillTable( "Table_Reconstruction", 3);
1040                 }
1041
1042                 TClonesArray* eNegPi0Dalitz = FindElectronFromPi0Dalitz(fEnegCandidateIndex, ::kElectron);
1043                 for(Int_t i=0; i < eNegPi0Dalitz->GetEntriesFast(); ++i)
1044                 {
1045                         AliKFParticle* eneg = (AliKFParticle*) eNegPi0Dalitz->At(i);
1046                         fHistograms->FillHistogram("MC_ESD_Pi0_EnegDalitz_Pt", eneg->GetPt());
1047                         fHistograms->FillHistogram("MC_ESD_Pi0_EnegDalitz_Eta", eneg->GetEta());
1048                         fHistograms->FillTable( "Table_Reconstruction", 4);
1049                 }
1050
1051                 TClonesArray* dalitzPairPi0 = FindDalitzPair(fEposCandidateIndex, fEnegCandidateIndex,1);
1052                 for(Int_t i=0; i < dalitzPairPi0->GetEntriesFast(); ++i)
1053                 {
1054                         TLorentzVector* dalitz = (TLorentzVector*) dalitzPairPi0->At(i);
1055                         fHistograms->FillHistogram("MC_ESD_Pi0_DalitzPair_Pt", dalitz->Pt());
1056                         fHistograms->FillHistogram("MC_ESD_Pi0_DalitzPair_Mass", dalitz->M());
1057                         fHistograms->FillHistogram( "Table_Reconstruction", 5 );
1058                 }
1059
1060                 TClonesArray* dalitzPairEta = FindDalitzPair(fEposCandidateIndex, fEnegCandidateIndex,2);
1061                 for(Int_t i=0; i < dalitzPairEta->GetEntriesFast(); ++i)
1062                 {
1063                         TLorentzVector* dalitz = (TLorentzVector*) dalitzPairEta->At(i);
1064                         fHistograms->FillHistogram("MC_ESD_Eta0_DalitzPair_Pt", dalitz->Pt());
1065                         fHistograms->FillHistogram("MC_ESD_Eta0_DalitzPair_InvMass", dalitz->M());
1066
1067                 }
1068
1069                 TClonesArray* lJpsiAll = FindJpsi(fEposCandidateIndex, fEnegCandidateIndex,-1);
1070
1071                 for(Int_t i=0; i < lJpsiAll->GetEntriesFast(); ++i)
1072                 {
1073                         TLorentzVector* jpsi = (TLorentzVector*) lJpsiAll->At(i);
1074                         fHistograms->FillHistogram("MC_ESD_Jpsi_Pt",jpsi->Pt());
1075                         fHistograms->FillHistogram("MC_ESD_Jpsi_InvMass",jpsi->M());
1076                         fHistograms->FillHistogram("MC_ESD_Jpsi_InvMass_vs_Pt",jpsi->M(),jpsi->Pt());
1077
1078                 }
1079
1080
1081                 TClonesArray* lJpsiChic0 = FindJpsi(fEposCandidateIndex, fEnegCandidateIndex,0);
1082
1083                 for(Int_t i=0; i < lJpsiChic0->GetEntriesFast(); ++i)
1084                 {
1085                         TLorentzVector* jpsi = (TLorentzVector*) lJpsiChic0->At(i);
1086                         fHistograms->FillHistogram("MC_ESD_Jpsi_Chic0_Pt",jpsi->Pt());
1087                         fHistograms->FillHistogram("MC_ESD_Jpsi_Chic0_InvMass",jpsi->M());
1088
1089                 }
1090                 TClonesArray* lJpsiChic1 = FindJpsi(fEposCandidateIndex, fEnegCandidateIndex,1);
1091
1092                 for(Int_t i=0; i < lJpsiChic1->GetEntriesFast(); ++i)
1093                 {
1094                         TLorentzVector* jpsi = (TLorentzVector*) lJpsiChic1->At(i);
1095                         fHistograms->FillHistogram("MC_ESD_Jpsi_Chic1_Pt",jpsi->Pt());
1096                         fHistograms->FillHistogram("MC_ESD_Jpsi_Chic1_InvMass",jpsi->M());
1097
1098                 }
1099                 TClonesArray* lJpsiChic2 = FindJpsi(fEposCandidateIndex, fEnegCandidateIndex,2);
1100
1101                 for(Int_t i=0; i < lJpsiChic2->GetEntriesFast(); ++i)
1102                 {
1103                         TLorentzVector* jpsi = (TLorentzVector*) lJpsiChic2->At(i);
1104                         fHistograms->FillHistogram("MC_ESD_Jpsi_Chic2_Pt",jpsi->Pt());
1105                         fHistograms->FillHistogram("MC_ESD_Jpsi_Chic2_InvMass",jpsi->M());
1106
1107                 }
1108
1109                 // psi pair for dalitz pairs
1110                 //if(fUsePsiPairCut)
1111                 FillPsiPair(ePosPi0Dalitz,eNegPi0Dalitz,"MC_ESD_Pi0_DalitzPair_PsiPair_vs_DPhi");
1112
1113                 delete ePosPi0Dalitz;
1114                 delete eNegPi0Dalitz;
1115                 delete dalitzPairPi0;
1116                 delete dalitzPairEta;
1117                 delete lJpsiAll;
1118                 delete lJpsiChic0;
1119                 delete lJpsiChic1;
1120                 delete lJpsiChic2;
1121                 // all gammas
1122                 TClonesArray* gamma = FindGamma(fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
1123                 for(Int_t i=0; i < gamma->GetEntriesFast(); ++i)
1124                 {
1125                         AliKFParticle* iGamma = (AliKFParticle*) gamma->At(i);
1126                         fHistograms->FillHistogram("MC_ESD_Gamma_Pt", iGamma->GetPt());
1127                         fHistograms->FillHistogram("MC_ESD_Gamma_Eta", iGamma->GetEta());
1128                 }
1129
1130                 delete gamma;
1131
1132                 // gamma from pi0 dalitz
1133                 TClonesArray* gammaPi0Dalitz = FindGammaFromPi0Dalitz(fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex);
1134                 for(Int_t i=0; i < gammaPi0Dalitz->GetEntriesFast(); ++i)
1135                 {
1136                         AliKFParticle* iGamma = (AliKFParticle*) gammaPi0Dalitz->At(i);
1137                         fHistograms->FillHistogram("MC_ESD_GammaPi0Dalitz_Pt", iGamma->GetPt());
1138                         fHistograms->FillHistogram("MC_ESD_GammaPi0Dalitz_Eta", iGamma->GetEta());
1139                         fHistograms->FillTable( "Table_Reconstruction", 6);
1140                 }
1141
1142                 delete gammaPi0Dalitz;
1143
1144                 TClonesArray* pi0Dalitz = FindParticleDalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,1);
1145
1146                 for(Int_t i=0; i < pi0Dalitz->GetEntriesFast(); ++i)
1147                 {
1148                         TLorentzVector* pi0 = (TLorentzVector*) pi0Dalitz->At(i);
1149
1150                         fHistograms->FillHistogram("MC_ESD_Pi0_P", pi0->P());
1151                         fHistograms->FillHistogram("MC_ESD_Pi0_Pt", pi0->Pt());
1152                         fHistograms->FillHistogram("MC_ESD_Pi0_Eta", pi0->Eta());
1153                         fHistograms->FillHistogram("MC_ESD_Pi0_Y", pi0->Rapidity());
1154                         fHistograms->FillHistogram("MC_ESD_Pi0_Phi", pi0->Phi());
1155                         fHistograms->FillHistogram("MC_ESD_Pi0_Y_vs_Pt",pi0->Pt(), pi0->Rapidity());
1156                         fHistograms->FillHistogram("MC_ESD_Pi0_InvMass", pi0->M());
1157                         fHistograms->FillHistogram("MC_ESD_Pi0_InvMass_vs_Pt",pi0->M(),pi0->Pt());
1158                         fHistograms->FillHistogram("MC_ESD_Pi0_InvMass_vs_Y", pi0->M(), pi0->Rapidity());
1159                         fHistograms->FillHistogram("MC_ESD_Pi0_InvMass_vs_Eta", pi0->M(),pi0->Eta());
1160                         fHistograms->FillHistogram( "Table_Reconstruction", 7);
1161                 }
1162                 delete pi0Dalitz;
1163
1164
1165                 TClonesArray* eta0Dalitz = FindParticleDalitz(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,2);
1166
1167
1168                 for(Int_t i=0; i < eta0Dalitz->GetEntriesFast(); ++i)
1169                 {
1170                         TLorentzVector* eta0 = (TLorentzVector*) eta0Dalitz->At(i);
1171
1172                         fHistograms->FillHistogram("MC_ESD_Eta0_P", eta0->P());
1173                         fHistograms->FillHistogram("MC_ESD_Eta0_Pt", eta0->Pt());
1174                         fHistograms->FillHistogram("MC_ESD_Eta0_Eta", eta0->Eta());
1175                         fHistograms->FillHistogram("MC_ESD_Eta0_Y", eta0->Rapidity());
1176                         fHistograms->FillHistogram("MC_ESD_Eta0_Phi", eta0->Phi());
1177                         fHistograms->FillHistogram("MC_ESD_Eta0_Pt_vs_Y", eta0->Pt(),eta0->Rapidity());
1178                         fHistograms->FillHistogram("MC_ESD_Eta0_InvMass", eta0->M());
1179                         fHistograms->FillHistogram("MC_ESD_Eta0_InvMass_vs_Pt", eta0->M(), eta0->Pt());
1180                         fHistograms->FillHistogram("MC_ESD_Eta0_InvMass_vs_Y", eta0->M(), eta0->Rapidity());
1181                         fHistograms->FillHistogram("MC_ESD_Eta0_InvMass_vs_Eta",eta0->M(),eta0->Eta());
1182                 }
1183                 delete eta0Dalitz;
1184
1185
1186                 TClonesArray* chic0Array = FindParticleChic(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,0);
1187
1188                 for(Int_t i=0; i < chic0Array->GetEntriesFast(); ++i)
1189                 {
1190                         TLorentzVector* chic0 = (TLorentzVector*) chic0Array->At(i);
1191                         fHistograms->FillHistogram("MC_ESD_Chic0_InvMass", chic0->M());
1192                         fHistograms->FillHistogram("MC_ESD_Chic0_InvMass_vs_Pt", chic0->M(),chic0->Pt());
1193                 }
1194                 delete chic0Array;
1195
1196                 TClonesArray* chic1Array = FindParticleChic(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,1);
1197
1198                 for(Int_t i=0; i < chic1Array->GetEntriesFast(); ++i)
1199                 {
1200                         TLorentzVector* chic1 = (TLorentzVector*) chic1Array->At(i);
1201                         fHistograms->FillHistogram("MC_ESD_Chic1_InvMass", chic1->M());
1202                         fHistograms->FillHistogram("MC_ESD_Chic1_InvMass_vs_Pt", chic1->M(), chic1->Pt());
1203                 }
1204                 delete chic1Array;
1205
1206                 TClonesArray* chic2Array = FindParticleChic(fEposCandidateIndex, fEnegCandidateIndex, fGammaCandidates, fGammaCandidatePosIndex, fGammaCandidateNegIndex,2);
1207
1208                 for(Int_t i=0; i < chic2Array->GetEntriesFast(); ++i)
1209                 {
1210                         TLorentzVector* chic2 = (TLorentzVector*) chic2Array->At(i);
1211                         fHistograms->FillHistogram("MC_ESD_Chic2_InvMass", chic2->M());
1212                         fHistograms->FillHistogram("MC_ESD_Chic2_InvMass_vs_Pt", chic2->M(), chic2->Pt());
1213                 }
1214                 delete chic2Array;
1215
1216
1217                 // psi pair for electrons from gamma conversions assuming they came from main vertex
1218                 // if(fUsePsiPairCut)
1219                 for(UInt_t i=0; i < fEposCandidateIndex.size(); ++i)
1220                 {
1221                         AliESDtrack* posTrack = fESDEvent->GetTrack(fEposCandidateIndex[i]);
1222                         Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1223
1224                         for(UInt_t j=0; j < fEnegCandidateIndex.size(); ++j)
1225                         {
1226                                 AliESDtrack* negTrack = fESDEvent->GetTrack(fEnegCandidateIndex[j]);
1227                                 Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1228
1229                                 if(!IsFromGammaConversion(posLabel,negLabel)) continue;
1230
1231                                 Double_t psiPair = GetPsiPair(posTrack, negTrack);
1232                                 Double_t deltaPhi = fMagFieldSign * TVector2::Phi_mpi_pi( negTrack->GetConstrainedParam()->Phi()-posTrack->GetConstrainedParam()->Phi());
1233
1234                                 fHistograms->FillHistogram("MC_ESD_EposEnegGamma_PsiPair_vs_DPhi", deltaPhi, psiPair);
1235                         }
1236                 }
1237                 // FIXME: eta -> e+e-gamma
1238         }
1239 }
1240
1241 //--------------------------------------------------------------------------
1242 Double_t AliAnalysisTaskGammaConvDalitz::Rapidity(const TParticle* p) const
1243 {
1244 //
1245 // Get rapidity
1246 //
1247         const double kEPSILON=1.e-16;
1248
1249         if(p->Energy() - TMath::Abs(p->Pz()) < kEPSILON )
1250         {
1251                 return 1.e10;
1252         }
1253         return 0.5*TMath::Log( (p->Energy()+p->Pz()) / (p->Energy()-p->Pz()) );
1254 }
1255
1256 //--------------------------------------------------------------------------
1257 void AliAnalysisTaskGammaConvDalitz::FillPsiPair(const TClonesArray* pos, const TClonesArray* neg, const TString& hName)
1258 {
1259 //
1260 // Fill histogram with psipair(pos,neg)
1261 //
1262         for(Int_t i=0; i < pos->GetEntriesFast(); ++i )
1263         {
1264                 AliKFParticle* posKF = (AliKFParticle*) pos->At(i);
1265                 for( Int_t j=0; j < neg->GetEntriesFast(); ++j )
1266                 {
1267                         AliKFParticle* negKF = (AliKFParticle*) neg->At(j);
1268                         Double_t psiPair = GetPsiPair(posKF, negKF);
1269                         Double_t deltaPhi = fMagFieldSign * TVector2::Phi_mpi_pi( negKF->GetPhi() - posKF->GetPhi());
1270                         fHistograms->FillHistogram(hName, deltaPhi, psiPair);
1271                 }
1272         }
1273 }
1274
1275 //--------------------------------------------------------------------------
1276 void AliAnalysisTaskGammaConvDalitz::FillAngle(const TClonesArray* x, const TClonesArray* y, const TString& hName)
1277 {
1278 //
1279 // Fill histogram with angle(x,y)
1280 //
1281         for(Int_t i=0; i < x->GetEntriesFast(); ++i )
1282         {
1283                 AliKFParticle* xKF = (AliKFParticle*) x->At(i);
1284                 TVector3 xMom(xKF->Px(),xKF->Py(),xKF->Pz());
1285                 for( Int_t j=0; j < y->GetEntriesFast(); ++j )
1286                 {
1287                         AliKFParticle* yKF = (AliKFParticle*) y->At(j);
1288                         TVector3 yMom(yKF->Px(),yKF->Py(),yKF->Pz());
1289                         fHistograms->FillHistogram(hName, xMom.Angle(yMom));
1290                 }
1291         }
1292 }
1293
1294 //--------------------------------------------------------------------------
1295 void AliAnalysisTaskGammaConvDalitz::FillPidTable(const TParticle* p, Int_t pid)
1296 {
1297 //
1298 // Fill table with pid info
1299 //
1300         Int_t iGen=-1;
1301         switch(TMath::Abs(p->GetPdgCode()))
1302         {
1303                 case ::kElectron:   iGen=0; break;
1304                 case ::kMuonMinus:  iGen=1; break;
1305                 case ::kPiPlus:     iGen=2; break;
1306                 case ::kKPlus:      iGen=3; break;
1307                 case ::kProton:     iGen=4; break;
1308                 default: iGen=-1;
1309         }
1310
1311         int jRec=-1;
1312         if(pid > -1 && pid < 5) jRec = pid;
1313
1314         if ((iGen > -1) && (jRec > -1))
1315         {
1316                 fHistograms->FillTable("Table_PID", iGen, jRec);
1317                 // sum
1318                 fHistograms->FillTable("Table_PID", iGen, 5);
1319                 fHistograms->FillTable("Table_PID", 5, jRec);
1320         }
1321 }
1322
1323 //--------------------------------------------------------------------------
1324 void AliAnalysisTaskGammaConvDalitz::GetGammaCandidates(TClonesArray*& gamma, vector<Int_t>& posIndex, vector<Int_t>& negIndex)
1325 {
1326 //
1327 // Make a copy of gamma candidates from V0reader
1328 //
1329         posIndex.clear();
1330         negIndex.clear();
1331
1332         if(gamma) delete gamma;
1333
1334         TClonesArray* gammaV0 = fV0Reader->GetCurrentEventGoodV0s();
1335
1336         gamma = new TClonesArray("AliKFParticle", gammaV0->GetEntriesFast());
1337         gamma->SetOwner(kTRUE);
1338
1339         // make a copy
1340         for(Int_t i=0; i < gammaV0->GetEntriesFast(); ++i)
1341         {
1342                 AliKFParticle* gamKF = (AliKFParticle*)gammaV0->At(i);
1343                 new ((*gamma)[i]) AliKFParticle(*gamKF);
1344                 posIndex.push_back(fV0Reader->GetPindex(i));
1345                 negIndex.push_back(fV0Reader->GetNindex(i));
1346         }
1347 }
1348
1349 //--------------------------------------------------------------------------
1350 TClonesArray* AliAnalysisTaskGammaConvDalitz::IndexToAliKFParticle(const vector<Int_t>& index, Int_t PDG)
1351 {
1352 //
1353 // Convert track index vector to AliKFParticle array
1354 //
1355         TClonesArray* indexKF = new TClonesArray("AliKFParticle",index.size());
1356         indexKF->SetOwner(kTRUE);
1357
1358         for(UInt_t i = 0; i < index.size(); ++i)
1359         {
1360                 AliESDtrack* t = fESDEvent->GetTrack(index[i]);
1361                 new((*indexKF)[i]) AliKFParticle(*t->GetConstrainedParam(), PDG);
1362         }
1363
1364         return indexKF;
1365 }
1366
1367 //--------------------------------------------------------------------------
1368 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindElectronFromPi0Dalitz(const vector<Int_t>& candidates, Int_t PDG)
1369 {
1370 //
1371 // Find true electrons from pi0 Dalitz decay candidates with MC
1372 //
1373         TClonesArray* elec = new TClonesArray("AliKFParticle");
1374         elec->SetOwner(kTRUE);
1375
1376         for(UInt_t i=0, j=0; i < candidates.size(); ++i)
1377         {
1378                 AliESDtrack* track = fESDEvent->GetTrack(candidates[i]);
1379                 Int_t trackLabel = TMath::Abs(track->GetLabel());
1380
1381                 if( fStack->Particle(trackLabel)->GetPdgCode() != PDG ) continue;
1382                 if( !IsPi0DalitzDaughter(trackLabel) ) continue;
1383
1384                 new ((*elec)[j++]) AliKFParticle(*track->GetConstrainedParam(), PDG);
1385         }
1386
1387         return elec;
1388 }
1389
1390 //--------------------------------------------------------------------------
1391 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindGammaFromPi0Dalitz(const TClonesArray* gamma, const vector<Int_t>& posIdx, const vector<Int_t>& negIdx)
1392 {
1393 //
1394 // Find true gammas from pi0 Dalitz decay candidates with MC
1395 //
1396         TClonesArray* gammaPi0 = new TClonesArray("AliKFParticle");
1397         gammaPi0->SetOwner(kTRUE);
1398
1399         for(Int_t i=0, j=0; i < gamma->GetEntriesFast(); ++i)
1400         {
1401                 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(i);
1402
1403                 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posIdx[i]))->GetLabel());
1404                 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negIdx[i]))->GetLabel());
1405
1406                 if( !HaveSameMother(labelv1,labelv2) ) continue;
1407
1408                 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
1409
1410                 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
1411
1412                 if( !IsPi0DalitzDaughter( labelGamma) ) continue;
1413
1414                 new ((*gammaPi0)[j++]) AliKFParticle(*gamKF);
1415         }
1416
1417         return gammaPi0;
1418 }
1419
1420 //--------------------------------------------------------------------------
1421 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindGamma(const TClonesArray* gamma, const vector<Int_t>& posIdx, const vector<Int_t>& negIdx)
1422 {
1423 //
1424 // Find true gammas from gamma candidates with MC
1425 //
1426         TClonesArray* gammaConv = new TClonesArray("AliKFParticle");
1427         gammaConv->SetOwner(kTRUE);
1428
1429         for(Int_t i=0, j=0; i < gamma->GetEntriesFast(); ++i)
1430         {
1431                 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(i);
1432
1433                 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posIdx[i]))->GetLabel());
1434                 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negIdx[i]))->GetLabel());
1435
1436                 if( !HaveSameMother(labelv1,labelv2) ) continue;
1437
1438                 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
1439
1440                 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
1441
1442                 new ((*gammaConv)[j++]) AliKFParticle(*gamKF);
1443         }
1444
1445         return gammaConv;
1446 }
1447
1448 //--------------------------------------------------------------
1449 void AliAnalysisTaskGammaConvDalitz::ESDtrackIndexCut(vector<Int_t>& pos, vector<Int_t>& neg, const TClonesArray* gamma)
1450 {
1451 //
1452 // Remove repeated electron candidate tracks
1453 // according to the gamma candidate array
1454 //
1455         vector<Bool_t> posTag(pos.size(),kTRUE);
1456         vector<Bool_t> negTag(neg.size(),kTRUE);
1457
1458         for(Int_t i=0; i < gamma->GetEntriesFast(); ++i)
1459         {
1460                 Int_t gamPosIndex = fGammaCandidatePosIndex[i];
1461                 Int_t gamNegIndex = fGammaCandidateNegIndex[i];
1462
1463                 for( UInt_t j=0; j < pos.size(); ++j )
1464                 {
1465                         if(pos[j] == gamPosIndex || pos[j] == gamNegIndex) posTag[j] = kFALSE;
1466                 }
1467
1468                 for( UInt_t j=0; j < neg.size(); ++j )
1469                 {
1470                         if(neg[j] == gamPosIndex || neg[j] == gamNegIndex) negTag[j] = kFALSE;
1471                 }
1472         }
1473
1474         CleanArray(pos, posTag);
1475         CleanArray(neg, negTag);
1476 }
1477
1478 //--------------------------------------------------------------------------
1479 void AliAnalysisTaskGammaConvDalitz::PsiPairCut(vector<Int_t>& pos, vector<Int_t>& neg)
1480 {
1481 //
1482 // Remove electron candidates from gamma conversions
1483 // according to the Psi pair angle
1484 //
1485     vector<Bool_t> posTag(pos.size(), kTRUE);
1486     vector<Bool_t> negTag(neg.size(), kTRUE);
1487
1488     for( UInt_t i=0; i < pos.size(); ++i )
1489     {
1490         AliESDtrack* posTrack = fESDEvent->GetTrack(pos[i]);
1491
1492         for( UInt_t j=0; j < neg.size(); ++j )
1493         {
1494             AliESDtrack* negTrack = fESDEvent->GetTrack(neg[j]);
1495
1496             Double_t psiPair = GetPsiPair(posTrack, negTrack);
1497             Double_t deltaPhi = fMagFieldSign * TVector2::Phi_mpi_pi( negTrack->GetConstrainedParam()->Phi()-posTrack->GetConstrainedParam()->Phi());
1498
1499             if(IsFromGammaConversion( psiPair, deltaPhi ))
1500             {
1501                 posTag[i] = kFALSE;
1502                 negTag[j] = kFALSE;
1503             }
1504         }
1505      }
1506
1507      CleanArray(pos, posTag);
1508      CleanArray(neg, negTag);
1509 }
1510
1511 //-----------------------------------------------------------------------------------
1512 void AliAnalysisTaskGammaConvDalitz::MassCut(vector<Int_t>& pos, vector<Int_t>& neg)
1513 {
1514 //
1515 // Remove electron candidates pairs
1516 // which have mass not in the range (fMassCutMin,fMassCutMax)
1517 //
1518     vector<Bool_t> posTag(pos.size(), kTRUE);
1519     vector<Bool_t> negTag(neg.size(), kTRUE);
1520
1521     for( UInt_t i=0; i < pos.size(); ++i )
1522     {
1523         AliESDtrack* posTrack = fESDEvent->GetTrack(pos[i]);
1524
1525         Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1526         TLorentzVector posLV;
1527         posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
1528
1529         for( UInt_t j=0; j < neg.size(); ++j )
1530         {
1531             AliESDtrack* negTrack = fESDEvent->GetTrack(neg[j]);
1532
1533             Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1534             TLorentzVector negLV;
1535             negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
1536
1537             TLorentzVector posnegLV = posLV + negLV;
1538
1539             if( (posnegLV.M() < fMassCutMin) || (posnegLV.M() > fMassCutMax) )
1540             {
1541                 posTag[i] = kFALSE;
1542                 negTag[j] = kFALSE;
1543             }
1544         }
1545      }
1546
1547      CleanArray(pos, posTag);
1548      CleanArray(neg, negTag);
1549 }
1550
1551 //-----------------------------------------------------------------------------------------------
1552 void AliAnalysisTaskGammaConvDalitz::CleanArray(vector<Int_t>& x, const vector<Bool_t>& tag)
1553 {
1554 //
1555 // Clean the x array according to the tag parameter
1556 //
1557         vector<Int_t> tmp;
1558
1559         for(UInt_t i=0; i< x.size(); ++i)
1560         {
1561                 if(tag[i]) tmp.push_back(x[i]);
1562         }
1563
1564         x = tmp;
1565 }
1566
1567 //--------------------------------------------------------------------------
1568 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)
1569 {
1570 //
1571 // Remove gamma candidates according to
1572 // the angle between the plane e+,e- and the gamma
1573 //
1574         vector<Bool_t> gammaTag(candidates->GetEntriesFast(), kTRUE);
1575
1576         for( UInt_t iPos=0; iPos < posIdx.size(); ++iPos )
1577         {
1578                 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1579                 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1580                 TVector3 xMom(posMom[0],posMom[1],posMom[2]);
1581
1582                 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1583                 {
1584                         AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1585                         Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1586                         TVector3 yMom(negMom[0],negMom[1],negMom[2]);
1587
1588                         // normal vector to x+y- plane
1589                         TVector3 planePosNeg = xMom.Cross(yMom);
1590                         for(Int_t i=0; i < candidates->GetEntriesFast(); ++i)
1591                         {
1592                                 AliKFParticle* gamKF = (AliKFParticle*)candidates->At(i);
1593                                 TVector3 gamMom(gamKF->Px(),gamKF->Py(),gamKF->Pz());
1594                                 if (planePosNeg.Angle(gamMom) < 1. || planePosNeg.Angle(gamMom) > 2.)
1595                                 {
1596                                         gammaTag[i] = kFALSE;
1597                                 }
1598                         }
1599                 }
1600         }
1601
1602         // Rebuild gamma candidates array
1603
1604         if(gamma) delete gamma;
1605         gamma = new TClonesArray("AliKFParticle");
1606         gamma->SetOwner(kTRUE);
1607
1608         posGamIdx.clear();
1609         negGamIdx.clear();
1610
1611         for(Int_t i=0, j=0; i < candidates->GetEntriesFast(); ++i)
1612         {
1613                 AliKFParticle* iGamma = (AliKFParticle*)candidates->At(i);
1614                 if(gammaTag[i])
1615                 {
1616                         new ((*gamma)[j++]) AliKFParticle(*iGamma);
1617                         posGamIdx.push_back(fV0Reader->GetPindex(i));
1618                         negGamIdx.push_back(fV0Reader->GetNindex(i));
1619                 }
1620         }
1621 }
1622
1623 //--------------------------------------------------------------------------
1624 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindDalitzPair(const TClonesArray* pos, const TClonesArray* neg)
1625 {
1626 //
1627 // Find Dalitz pair candidates
1628 //
1629         TClonesArray* dalitz = new TClonesArray("TLorentzVector");
1630         dalitz->SetOwner(kTRUE);
1631
1632         for( Int_t iPos=0, j=0; iPos < pos->GetEntriesFast(); ++iPos )
1633         {
1634                 AliKFParticle* posKF = (AliKFParticle*)pos->At(iPos);
1635
1636                 TLorentzVector posLV;
1637                 posLV.SetXYZM(posKF->Px(),posKF->Py(),posKF->Pz(),fkElectronMass);
1638
1639                 for( Int_t iNeg=0; iNeg < neg->GetEntriesFast(); ++iNeg )
1640                 {
1641                         AliKFParticle* negKF = (AliKFParticle*)neg->At(iNeg);
1642
1643                         TLorentzVector negLV;
1644                         negLV.SetXYZM(negKF->Px(),negKF->Py(),negKF->Pz(),fkElectronMass);
1645
1646                         if(fUseAliKF)
1647                         {
1648                                 AliKFParticle posNeg( *posKF, *negKF);
1649
1650                                 TLorentzVector posNegLV;
1651                                 posNegLV.SetXYZM(posNeg.Px(), posNeg.Py(), posNeg.Pz(), posNeg.GetMass());
1652                                 new ((*dalitz)[j++]) TLorentzVector(posNegLV);
1653                         }
1654                         else
1655                         {
1656                                 new ((*dalitz)[j++]) TLorentzVector(posLV + negLV);
1657                         }
1658                 }
1659         }
1660
1661         return dalitz;
1662 }
1663
1664 //--------------------------------------------------------------------------
1665 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindParticleDalitz(const TClonesArray* pos, const TClonesArray* neg, const TClonesArray* gamma,Int_t opc)
1666 {
1667 //
1668 // Find pi0 Dalitz decay candidates
1669 //
1670         TClonesArray* pi0 = new TClonesArray("TLorentzVector");
1671         pi0->SetOwner(kTRUE);
1672
1673         for( Int_t iPos=0, j=0; iPos < pos->GetEntriesFast(); ++iPos )
1674         {
1675                 AliKFParticle* posKF = (AliKFParticle*)pos->At(iPos);
1676
1677                 TLorentzVector posLV;
1678                 posLV.SetXYZM(posKF->Px(),posKF->Py(),posKF->Pz(),fkElectronMass);
1679
1680                 for( Int_t iNeg=0; iNeg < neg->GetEntriesFast(); ++iNeg )
1681                 {
1682                         AliKFParticle* negKF = (AliKFParticle*)neg->At(iNeg);
1683
1684                         TLorentzVector negLV;
1685                         negLV.SetXYZM(negKF->Px(),negKF->Py(),negKF->Pz(),fkElectronMass);
1686
1687                         AliKFParticle posNegKF(*posKF,*negKF);
1688
1689
1690                         for(Int_t iGam=0; iGam < gamma->GetEntriesFast(); ++iGam)
1691                         {
1692                                 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(iGam);
1693                                 AliKFParticle posNegGam( *posKF, *negKF, *gamKF );
1694
1695                                 Double_t lDiffMass = posNegGam.GetMass() - posNegKF.GetMass();
1696
1697                                 if( opc == 1 )
1698                                 {
1699                                     fHistograms->FillHistogram("ESD_BKG_PrevGamma_InvMass_Diff",lDiffMass );
1700                                     fHistograms->FillHistogram("ESD_BKG_PrevGamma_InvMass_vs_Pt_Diff",lDiffMass,posNegGam.GetPt());
1701                                 }
1702                                 else if ( opc == 2 )
1703                                 {
1704                                     fHistograms->FillHistogram("ESD_BKG_BGHandler_InvMass_Diff",lDiffMass );
1705                                     fHistograms->FillHistogram("ESD_BKG_BGHandler_InvMass_vs_Pt_Diff",lDiffMass,posNegGam.GetPt());
1706                                 }
1707                                 else if ( opc == 3 )
1708                                 {
1709                                     fHistograms->FillHistogram("ESD_BKG_Electron_InvMass_Diff",lDiffMass );
1710                                     fHistograms->FillHistogram("ESD_BKG_Electron_InvMass_vs_Pt_Diff",lDiffMass,posNegGam.GetPt());
1711                                 }
1712
1713                                 if(fUseAliKF)
1714                                 {
1715
1716                                         TLorentzVector posNegGamLV;
1717                                         posNegGamLV.SetXYZM(posNegGam.Px(),posNegGam.Py(),posNegGam.Pz(),posNegGam.GetMass());
1718                                         new ((*pi0)[j++]) TLorentzVector(posNegGamLV);
1719                                 }
1720                                 else
1721                                 {
1722                                         TLorentzVector gamLV;
1723                                         gamLV.SetXYZM(gamKF->Px(),gamKF->Py(),gamKF->Pz(),0);
1724
1725                                         new ((*pi0)[j++]) TLorentzVector(posLV + negLV + gamLV);
1726                                 }
1727                         }
1728                 }
1729         }
1730
1731         return pi0;
1732 }
1733
1734 //--------------------------------------------------------------------------
1735 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindDalitzPair(const vector<Int_t>& posIdx, const vector<Int_t>& negIdx,Int_t motherOpc)
1736 {
1737 //
1738 // Find true Dalitz pairs from Dalitz pair candidats with MC
1739 //
1740         TClonesArray* dalitz = new TClonesArray("TLorentzVector");
1741         dalitz->SetOwner(kTRUE);
1742
1743         for( UInt_t iPos=0, j=0; iPos < posIdx.size(); ++iPos )
1744         {
1745                 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1746                 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1747                 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1748
1749                 TLorentzVector posLV;
1750                 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
1751
1752                 AliKFParticle posKF( *posTrack->GetConstrainedParam(), ::kPositron );
1753
1754                 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1755                 {
1756                         AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1757                         Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1758
1759                         if(!IsDalitzPair(posLabel,negLabel,motherOpc)) continue;
1760
1761                         if(fUseAliKF)
1762                         {
1763                                 AliKFParticle negKF( *negTrack->GetConstrainedParam(), ::kElectron );
1764                                 AliKFParticle posNeg( posKF, negKF);
1765
1766                                 TLorentzVector posNegLV;
1767                                 posNegLV.SetXYZM(posNeg.Px(),posNeg.Py(),posNeg.Pz(),posNeg.GetMass());
1768
1769                                 new ((*dalitz)[j++]) TLorentzVector(posNegLV);
1770                         }
1771                         else // TLorentzVector
1772                         {
1773                                 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1774
1775                                 TLorentzVector negLV;
1776                                 negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
1777
1778                                 new ((*dalitz)[j++]) TLorentzVector(posLV + negLV);
1779                         }
1780                 }
1781         }
1782
1783         return dalitz;
1784 }
1785
1786 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindJpsi(const vector<Int_t>& posIdx, const vector<Int_t>& negIdx,Int_t motherOpc)
1787 {
1788 //
1789 // Find true Jpsi's
1790 // If mother
1791 // -1: from the all sources
1792 // 0: from the Chic_0
1793 // 1: from the Chic_1
1794 // 2: from the Chic_2
1795
1796         TClonesArray* jPsi = new TClonesArray("TLorentzVector");
1797         jPsi->SetOwner(kTRUE);
1798
1799         for( UInt_t iPos=0, j=0; iPos < posIdx.size(); ++iPos )
1800         {
1801                 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1802                 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1803                 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1804
1805                 if( fStack->Particle(posLabel)->GetPdgCode() != ::kPositron ) continue;
1806
1807                 TLorentzVector posLV;
1808                 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
1809
1810                 AliKFParticle posKF( *posTrack->GetConstrainedParam(), ::kPositron );
1811
1812                 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1813                 {
1814                         AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1815                         Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1816
1817                         if( fStack->Particle(negLabel)->GetPdgCode() != ::kElectron ) continue;
1818
1819                         if( !HaveSameMother(posLabel,negLabel) ) continue;
1820
1821                         Int_t motherLabel = fStack->Particle(negLabel)->GetMother(0);
1822                         TParticle *motherJpsi = fStack->Particle(motherLabel);
1823
1824                         if( motherJpsi->GetPdgCode() != 443 ){
1825                             continue;
1826                         }
1827
1828
1829
1830                         if( motherOpc > -1 )
1831                         {
1832                              if( motherJpsi->GetMother(0) < 0 ) continue;
1833
1834                              TParticle *gmotherChic = fStack->Particle(motherJpsi->GetMother(0));
1835                              Int_t pdgCode = gmotherChic->GetPdgCode();
1836
1837                              Bool_t lson = kTRUE;
1838
1839                              switch(motherOpc){
1840
1841                                     case 0:     if ( pdgCode != 10441 )
1842                                                 lson = kFALSE;
1843                                                 break;
1844                                     case 1:     if ( pdgCode != 20443 )
1845                                                 lson = kFALSE;
1846                                                 break;
1847                                     case 2:     if ( pdgCode != 445   )
1848                                                 lson = kFALSE;
1849                                                 break;
1850                              }
1851
1852                             if( lson == kFALSE ) continue;
1853                         }
1854
1855
1856                         if(fUseAliKF)
1857                         {
1858                                 AliKFParticle negKF( *negTrack->GetConstrainedParam(), ::kElectron );
1859                                 AliKFParticle posNeg( posKF, negKF);
1860
1861                                 TLorentzVector posNegLV;
1862                                 posNegLV.SetXYZM(posNeg.Px(),posNeg.Py(),posNeg.Pz(),posNeg.GetMass());
1863
1864                                 new ((*jPsi)[j++]) TLorentzVector(posNegLV);
1865                         }
1866                         else // TLorentzVector
1867                         {
1868                                 Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1869
1870                                 TLorentzVector negLV;
1871                                 negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
1872
1873                                 new ((*jPsi)[j++]) TLorentzVector(posLV + negLV);
1874                         }
1875                 }
1876         }
1877
1878         return jPsi;
1879 }
1880
1881
1882
1883 //--------------------------------------------------------------------------
1884 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindParticleDalitz(const vector<Int_t>& posIdx, const vector<Int_t>& negIdx, const TClonesArray* gamma, const vector<Int_t>& posGam, const vector<Int_t>& negGam,Int_t motherOpc)
1885 {
1886 //
1887 // Find true pi0 Dalitz decay from pi0 candidates with MC
1888 //
1889         TClonesArray* pi0 = new TClonesArray("TLorentzVector");
1890         pi0->SetOwner(kTRUE);
1891
1892         for( UInt_t iPos=0, j=0; iPos < posIdx.size(); ++iPos )
1893         {
1894                 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1895                 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1896                 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1897
1898                 TLorentzVector posLV;
1899                 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
1900
1901                 AliKFParticle posKF( *posTrack->GetConstrainedParam(), ::kPositron );
1902
1903                 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1904                 {
1905                         AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1906                         Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1907
1908
1909                         if( !HaveSameMother(posLabel,negLabel) ) continue; //Check if both particles have same mother
1910                         if(!IsDalitzPair(posLabel,negLabel,motherOpc)) continue; //check if mohter is eta0 or pi0
1911
1912                         Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
1913
1914                         TLorentzVector negLV;
1915                         negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
1916
1917                         AliKFParticle negKF( *negTrack->GetConstrainedParam(), ::kElectron );
1918
1919                         for(Int_t iGam=0; iGam < gamma->GetEntriesFast(); ++iGam)
1920                         {
1921                                 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(iGam);
1922
1923                                 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posGam[iGam]))->GetLabel());
1924                                 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negGam[iGam]))->GetLabel());
1925
1926                                 if( !HaveSameMother(labelv1,labelv2) ) continue;
1927
1928                                 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
1929
1930                                 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
1931
1932
1933                                 if( !HaveSameMother(labelGamma, posLabel) ) continue;
1934
1935
1936                                 if(fUseAliKF)
1937                                 {
1938                                         AliKFParticle posNegGam( posKF, negKF, *gamKF );
1939                                         TLorentzVector posNegGamLV;
1940                                         posNegGamLV.SetXYZM(posNegGam.Px(),posNegGam.Py(),posNegGam.Pz(),posNegGam.GetMass());
1941                                         new ((*pi0)[j++]) TLorentzVector(posNegGamLV);
1942                                 }
1943                                 else // TLorentzVector
1944                                 {
1945                                         TLorentzVector gamLV;
1946                                         gamLV.SetXYZM(gamKF->Px(),gamKF->Py(),gamKF->Pz(),0);
1947
1948                                         new ((*pi0)[j++]) TLorentzVector(posLV + negLV + gamLV);
1949                                 }
1950                         }
1951                 }
1952         }
1953
1954         return pi0;
1955 }
1956 TClonesArray* AliAnalysisTaskGammaConvDalitz::FindParticleChic(const vector<Int_t>& posIdx, const vector<Int_t>& negIdx, const TClonesArray* gamma, const vector<Int_t>& posGam, const vector<Int_t>& negGam,Int_t motherOpc)
1957 {
1958 //
1959 // Find true pi0 Dalitz decay from pi0 candidates with MC
1960 //
1961         TClonesArray* chic = new TClonesArray("TLorentzVector");
1962         chic->SetOwner(kTRUE);
1963
1964         for( UInt_t iPos=0, j=0; iPos < posIdx.size(); ++iPos )
1965         {
1966                 AliESDtrack* posTrack = fESDEvent->GetTrack(posIdx[iPos]);
1967                 Double_t posMom[3]; posTrack->GetConstrainedPxPyPz(posMom);
1968                 Int_t posLabel = TMath::Abs(posTrack->GetLabel());
1969
1970                 if( fStack->Particle(posLabel)->GetPdgCode() != ::kPositron ) continue;
1971
1972                 TLorentzVector posLV;
1973                 posLV.SetXYZM(posMom[0],posMom[1],posMom[2],fkElectronMass);
1974
1975                 AliKFParticle posKF( *posTrack->GetConstrainedParam(), ::kPositron );
1976
1977                 for( UInt_t iNeg=0; iNeg < negIdx.size(); ++iNeg )
1978                 {
1979                         AliESDtrack* negTrack = fESDEvent->GetTrack(negIdx[iNeg]);
1980                         Int_t negLabel = TMath::Abs(negTrack->GetLabel());
1981
1982
1983                         if( fStack->Particle(negLabel)->GetPdgCode() != ::kElectron ) continue;
1984
1985
1986                         if( !HaveSameMother(posLabel,negLabel) ) continue;
1987
1988                         Int_t jpsiLabel = fStack->Particle(negLabel)->GetMother(0);
1989
1990                         if( fStack->Particle(jpsiLabel)->GetPdgCode() != 443 ) continue;
1991
1992                         TParticle *jpsiParticle = fStack->Particle(jpsiLabel);
1993
1994                         if ( jpsiParticle->GetMother(0) < 0 ) continue;
1995
1996
1997                         Int_t chicLabel = jpsiParticle->GetMother(0);
1998
1999                         Int_t pdgCode = fStack->Particle(chicLabel)->GetPdgCode();
2000
2001
2002                         Bool_t lSon = kTRUE;
2003
2004                         switch(motherOpc){
2005
2006                                     case 0:     if ( pdgCode != 10441 )
2007                                                 lSon = kFALSE;
2008                                                 break;
2009                                     case 1:     if ( pdgCode != 20443 )
2010                                                 lSon = kFALSE;
2011                                                 break;
2012                                     case 2:     if ( pdgCode != 445   )
2013                                                 lSon = kFALSE;
2014                                                 break;
2015                         }
2016
2017
2018                        if( lSon == kFALSE ) continue;
2019
2020
2021
2022
2023
2024                         Double_t negMom[3]; negTrack->GetConstrainedPxPyPz(negMom);
2025
2026                         TLorentzVector negLV;
2027                         negLV.SetXYZM(negMom[0],negMom[1],negMom[2],fkElectronMass);
2028
2029                         AliKFParticle negKF( *negTrack->GetConstrainedParam(), ::kElectron );
2030
2031                         for(Int_t iGam=0; iGam < gamma->GetEntriesFast(); ++iGam)
2032                         {
2033                                 AliKFParticle* gamKF = (AliKFParticle*)gamma->At(iGam);
2034
2035                                 Int_t labelv1 = TMath::Abs((fESDEvent->GetTrack(posGam[iGam]))->GetLabel());
2036                                 Int_t labelv2 = TMath::Abs((fESDEvent->GetTrack(negGam[iGam]))->GetLabel());
2037
2038                                 if( !HaveSameMother(labelv1,labelv2) ) continue;
2039
2040                                 Int_t labelGamma = TMath::Abs(fStack->Particle(labelv1)->GetMother(0));
2041
2042                                 if( fStack->Particle(labelGamma)->GetPdgCode() != ::kGamma ) continue;
2043
2044
2045                                 if( !HaveSameMother(labelGamma, jpsiLabel) ) continue;
2046
2047
2048                                 if(fUseAliKF)
2049                                 {
2050                                         AliKFParticle posNegGam( posKF, negKF, *gamKF );
2051                                         TLorentzVector posNegGamLV;
2052                                         posNegGamLV.SetXYZM(posNegGam.Px(),posNegGam.Py(),posNegGam.Pz(),posNegGam.GetMass());
2053                                         new ((*chic)[j++]) TLorentzVector(posNegGamLV);
2054                                 }
2055                                 else // TLorentzVector
2056                                 {
2057                                         TLorentzVector gamLV;
2058                                         gamLV.SetXYZM(gamKF->Px(),gamKF->Py(),gamKF->Pz(),0);
2059
2060                                         new ((*chic)[j++]) TLorentzVector(posLV + negLV + gamLV);
2061                                 }
2062                         }
2063                 }
2064         }
2065
2066         return chic;
2067 }
2068
2069 //-----------------------------------------------------------------------------------------------
2070 void AliAnalysisTaskGammaConvDalitz::UpdateGammaPool(const TClonesArray* gamma)
2071 {
2072 //
2073 // Update gamma event pool for background computation
2074 //
2075         if( fDebug ) AliInfo("=> UpdateGammaPool");
2076
2077         // cycle
2078         for(Int_t j=0; j< gamma->GetEntriesFast(); ++j)
2079         {
2080                 if((AliKFParticle*)fGammaPool->At(fGamPoolPos)) delete (AliKFParticle*)fGammaPool->RemoveAt(fGamPoolPos);
2081                 new ((*fGammaPool)[fGamPoolPos]) AliKFParticle( *((AliKFParticle*)gamma->At(j)));
2082                 ++fGamPoolPos;
2083                 if(fGamPoolPos == fPoolMaxSize)
2084                 {
2085                         fGamPoolPos = 0;
2086                 }
2087         }
2088 }
2089
2090 void AliAnalysisTaskGammaConvDalitz::UpdateElectronPool(TClonesArray* elec) // FIXME: const
2091 {
2092 //
2093 // Update electron event pool for background computation
2094 //
2095         Int_t multiplicity = fV0Reader->CountESDTracks();
2096
2097
2098         fBGEventHandler->AddElectronEvent(elec,fESDEvent->GetPrimaryVertex()->GetZ(),multiplicity);
2099 }
2100
2101 //-----------------------------------------------------------------------------------------------
2102 TClonesArray* AliAnalysisTaskGammaConvDalitz::GammasFromBGHandler() const
2103 {
2104 //
2105 // Gamma copy from events with same multiplicity and Z
2106 //
2107         if( fDebug ) AliInfo("=> GammasFromBGHandler");
2108
2109         Int_t zbin = fBGEventHandler->GetZBinIndex(fV0Reader->GetVertexZ());
2110         Int_t mbin = fBGEventHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2111
2112
2113
2114         TClonesArray* gammaPool = new TClonesArray("AliKFParticle");
2115         gammaPool->SetOwner(kTRUE);
2116
2117         for( Int_t iEventBG=0; iEventBG < fV0Reader->GetNBGEvents(); ++iEventBG )
2118         {
2119                 AliGammaConversionKFVector* gammaV0s = fBGEventHandler->GetBGGoodV0s(zbin,mbin,iEventBG);
2120                 for( UInt_t i = 0; i < gammaV0s->size(); ++i)
2121                 {
2122                         new ((*gammaPool)[i]) AliKFParticle( *((AliKFParticle*)gammaV0s->at(i)) );
2123                 }
2124         }
2125
2126         return gammaPool;
2127 }
2128
2129 //-----------------------------------------------------------------------------------------------
2130 TClonesArray* AliAnalysisTaskGammaConvDalitz::ElectronFromBGHandler() const
2131 {
2132 //
2133 // Electron copy from events with same multiplicity and Z
2134 //
2135         if( fDebug ) AliInfo("=> ElectronFromBGHandler");
2136
2137         TClonesArray* electronPool = new TClonesArray("AliKFParticle");
2138         electronPool->SetOwner(kTRUE);
2139
2140         Int_t multiplicity = fV0Reader->CountESDTracks();
2141
2142
2143
2144
2145         for( Int_t iEventBG=0; iEventBG < fV0Reader->GetNBGEvents(); ++iEventBG )
2146         {
2147                 AliGammaConversionKFVector* electronNeg =  fBGEventHandler->GetBGGoodENeg(iEventBG,fESDEvent->GetPrimaryVertex()->GetZ(),multiplicity);
2148                 for (UInt_t i = 0; i < electronNeg->size(); ++i  )
2149                 {
2150                         new ((*electronPool)[i]) AliKFParticle( *((AliKFParticle*)electronNeg->at(i)) );
2151                 }
2152         }
2153
2154         return electronPool;
2155 }
2156
2157 //-----------------------------------------------------------------------------------------------
2158 Int_t AliAnalysisTaskGammaConvDalitz::GetMonteCarloPid(const AliESDtrack* t) const
2159 {
2160 //
2161 // Get track pid according to MC
2162 //
2163     Int_t label   = TMath::Abs(t->GetLabel());
2164     Int_t pdgCode = TMath::Abs(fStack->Particle(label)->GetPdgCode());
2165
2166     switch(pdgCode)
2167     {
2168          case ::kElectron:  return AliPID::kElectron;
2169          case ::kMuonMinus: return AliPID::kMuon;
2170          case ::kPiPlus:    return AliPID::kPion;
2171          case ::kKPlus:     return AliPID::kKaon;
2172          case ::kProton:    return AliPID::kProton;
2173     }
2174
2175     return -1;
2176 }
2177
2178 //-----------------------------------------------------------------------------------------------
2179 //FIXME PID ITS
2180 // NOTE prior should be estimated from data
2181 // NOTE: move to config
2182
2183 Int_t AliAnalysisTaskGammaConvDalitz::GetBayesPid(const AliESDtrack* t, Int_t trackType ) const
2184 {
2185 //
2186 // Get track pid according to Bayes' formula
2187 //
2188         double priors[AliPID::kSPECIES] = {0.009, 0.01, 0.82, 0.10, 0.05};
2189         Double_t detectoProb[AliPID::kSPECIES];
2190
2191         if( trackType == kITSsaTrack )  // ITS standalone pid
2192         {
2193            t->GetITSpid( detectoProb );
2194         }
2195         else  // global track
2196         {
2197            t->GetESDpid( detectoProb );
2198         }
2199
2200         AliPID bayesPID( detectoProb );
2201         return bayesPID.GetMostProbable( priors );
2202 }
2203
2204 //-----------------------------------------------------------------------------------------------
2205 Int_t AliAnalysisTaskGammaConvDalitz::GetNSigmaPid(const AliESDtrack* t, Int_t trackType ) const
2206 {
2207 //
2208 // Get track pid according to a n-sigma cut around ITS and/or TPC signals
2209 //
2210         if( trackType == kITSsaTrack)   // ITS standalone tracks
2211         {
2212             Double_t mom = t->GetP();
2213
2214
2215             // ITS Number of sigmas (FIXME: add new fESDpidCuts)
2216             // NOTE there is not AliESDpidCuts::SetITSnSigmaCut yet
2217             Double_t nElecSigma = fESDpid->NumberOfSigmasITS(t, AliPID::kElectron );
2218             Double_t nPionSigma = fESDpid->NumberOfSigmasITS(t, AliPID::kPion );
2219
2220             if( nElecSigma < 4. && nElecSigma > -3. && mom < .2  &&  nPionSigma < -1.5 )
2221             {
2222                 return AliPID::kElectron;
2223             }
2224         }
2225         else // global track
2226         {
2227             Double_t nElecSigma   = fESDpid->NumberOfSigmasTPC(t, AliPID::kElectron );
2228             Double_t nPionSigma   = fESDpid->NumberOfSigmasTPC(t, AliPID::kPion );
2229             Double_t nKaonSigma   = TMath::Abs(fESDpid->NumberOfSigmasTPC(t, AliPID::kKaon ));
2230             Double_t nProtonSigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(t, AliPID::kProton));
2231             if(     nElecSigma   > fNSigmaBelowElecTPCbethe && nElecSigma < fNSigmaAboveElecTPCbethe &&
2232                     nPionSigma   > fNSigmaAbovePionTPCbethe && //NOTE mom > 0.5
2233                     nKaonSigma   > fNSigmaAboveKaonTPCbethe &&
2234                     nProtonSigma > fNSigmaAboveProtonTPCbethe )
2235             {
2236                return AliPID::kElectron;
2237             }
2238             // NOTE: add other particle types
2239           }
2240
2241           return -1;
2242 }
2243
2244 //-----------------------------------------------------------------------------------------------
2245 Bool_t AliAnalysisTaskGammaConvDalitz::IsDalitzPair( Int_t posLabel, Int_t negLabel,Int_t motherOpc ) const
2246 {
2247 //
2248 // Returns true if the two particles is a Dalitz pair
2249 //
2250 //motherOpc 1:  for Pi0Dalitz
2251 //motherOpc 2:  for EtaDalitz
2252
2253         if(!HaveSameMother(posLabel, negLabel)) return kFALSE;
2254
2255         TParticle* pos = fStack->Particle( posLabel );
2256         TParticle* neg = fStack->Particle( negLabel );
2257
2258         if( pos->GetPdgCode() != ::kPositron ) return kFALSE;
2259         if( neg->GetPdgCode() != ::kElectron ) return kFALSE;
2260
2261         //if( pos->GetUniqueID() != ::kPDecay ) return kFALSE;
2262         //if( neg->GetUniqueID() != ::kPDecay ) return kFALSE;
2263
2264         Int_t motherLabel = pos->GetMother(0);
2265         if( motherLabel < 0 ) return kFALSE;
2266
2267         TParticle* mother = fStack->Particle( motherLabel );
2268
2269         if( mother->GetNDaughters() != 3)    return kFALSE;
2270
2271         if( motherOpc == 1 ){ //Pi0Dalitz
2272                if( mother->GetPdgCode() != ::kPi0 ) return kFALSE;
2273         }
2274         else if(motherOpc == 2){
2275                if( mother->GetPdgCode() != 221 ) return kFALSE;
2276         }
2277            else {
2278                 return kFALSE;
2279         }
2280         // NOTE: one of them must be a photon
2281
2282         return kTRUE;
2283 }
2284
2285 //-----------------------------------------------------------------------------------------------
2286 Bool_t AliAnalysisTaskGammaConvDalitz::IsPi0DalitzDaughter( Int_t label ) const
2287 {
2288 //
2289 // Returns true if the particle comes from Pi0 -> e+ e- gamma
2290 //
2291         Bool_t ePlusFlag  = kFALSE;
2292         Bool_t eMinusFlag = kFALSE;
2293         Bool_t gammaFlag  = kFALSE;
2294
2295         Int_t motherLabel = fStack->Particle( label )->GetMother(0);
2296
2297         if( motherLabel < 0 ) return kFALSE;
2298
2299         TParticle* mother = fStack->Particle( motherLabel );
2300
2301         if ( mother->GetPdgCode() != ::kPi0 ) return kFALSE;
2302
2303         if ( mother->GetNDaughters() != 3 ) return kFALSE;
2304
2305         for( Int_t idx = mother->GetFirstDaughter(); idx <= mother->GetLastDaughter(); ++idx )
2306         {
2307                 switch( fStack->Particle(idx)->GetPdgCode())
2308                 {
2309                         case ::kPositron:
2310                                 ePlusFlag  = kTRUE;
2311                                 break;
2312                         case ::kElectron:
2313                                 eMinusFlag = kTRUE;
2314                                 break;
2315                         case ::kGamma:
2316                                 gammaFlag  = kTRUE;
2317                                 break;
2318                 }
2319         }
2320
2321         return ( ePlusFlag && eMinusFlag && gammaFlag );
2322 }
2323
2324 //--------------------------------------------------------------------------
2325 Bool_t AliAnalysisTaskGammaConvDalitz::IsFromGammaConversion( Double_t psiPair, Double_t deltaPhi ) const
2326 {
2327 //
2328 // Returns true if it is a gamma conversion according to psi pair value
2329 //
2330         return ( (deltaPhi > fDeltaPhiCutMin  &&  deltaPhi < fDeltaPhiCutMax) &&
2331         TMath::Abs(psiPair) < ( fPsiPairCut - fPsiPairCut/fDeltaPhiCutMax * deltaPhi ) );
2332 }
2333
2334 //--------------------------------------------------------------------------
2335 Bool_t AliAnalysisTaskGammaConvDalitz::IsFromGammaConversion( Int_t posLabel, Int_t negLabel ) const
2336 {
2337 //
2338 // Returns true if it is a gamma conversion according to MC
2339 //
2340         if( !HaveSameMother(posLabel,negLabel) ) return kFALSE;
2341
2342         TParticle* pos = fStack->Particle( posLabel );
2343         TParticle* neg = fStack->Particle( negLabel );
2344
2345         if( pos->GetPdgCode() != ::kPositron ) return kFALSE;
2346         if( neg->GetPdgCode() != ::kElectron ) return kFALSE;
2347
2348         if( pos->GetUniqueID() != kPPair ) return kFALSE;
2349
2350         Int_t motherLabel = pos->GetMother(0);
2351         if( motherLabel < 0 ) return kFALSE;
2352
2353         TParticle* mother = fStack->Particle( motherLabel );
2354
2355         if( mother->GetPdgCode() != ::kGamma ) return kFALSE;
2356
2357         return kTRUE;
2358 }
2359
2360 //-----------------------------------------------------------------------------------------------
2361 Bool_t AliAnalysisTaskGammaConvDalitz::HaveSameMother( Int_t label1, Int_t label2 ) const
2362 {
2363 //
2364 // Returns true if the two particle have the same mother
2365 //
2366         if(fStack->Particle( label1 )->GetMother(0) < 0 ) return kFALSE;
2367         return (fStack->Particle( label1 )->GetMother(0) == fStack->Particle( label2 )->GetMother(0));
2368 }
2369
2370 //-----------------------------------------------------------------------------------------------
2371 Double_t AliAnalysisTaskGammaConvDalitz::GetPsiPair( const AliESDtrack* trackPos, const AliESDtrack* trackNeg ) const
2372 {
2373 //
2374 // This angle is a measure for the contribution of the opening in polar
2375 // direction Δ0 to the opening angle ξ Pair
2376 //
2377 // Ref. Measurement of photons via conversion pairs with the PHENIX experiment at RHIC
2378 //      Master Thesis. Thorsten Dahms. 2005
2379 // https://twiki.cern.ch/twiki/pub/ALICE/GammaPhysicsPublications/tdahms_thesis.pdf
2380 //
2381     Double_t momPos[3];
2382     Double_t momNeg[3];
2383     if( trackPos->GetConstrainedPxPyPz(momPos) == 0 ) trackPos->GetPxPyPz( momPos );
2384     if( trackNeg->GetConstrainedPxPyPz(momNeg) == 0 ) trackNeg->GetPxPyPz( momNeg );
2385
2386     TVector3 posDaughter;
2387     TVector3 negDaughter;
2388
2389     posDaughter.SetXYZ( momPos[0], momPos[1], momPos[2] );
2390     negDaughter.SetXYZ( momNeg[0], momNeg[1], momNeg[2] );
2391
2392     Double_t deltaTheta = negDaughter.Theta() - posDaughter.Theta();
2393     Double_t openingAngle =  posDaughter.Angle( negDaughter );  //TMath::ACos( posDaughter.Dot(negDaughter)/(negDaughter.Mag()*posDaughter.Mag()) );
2394     if( openingAngle < 1e-20 ) return 0.;
2395     Double_t psiAngle = TMath::ASin( deltaTheta/openingAngle );
2396
2397     return psiAngle;
2398 }
2399
2400 //-----------------------------------------------------------------------------------------------
2401 Double_t AliAnalysisTaskGammaConvDalitz::GetPsiPair(const AliKFParticle* xPos, const AliKFParticle* yNeg ) const
2402 {
2403 //
2404 // Get psi pair value
2405 //
2406     TVector3 pos(xPos->GetPx(), xPos->GetPy(), xPos->GetPz());
2407     TVector3 neg(yNeg->GetPx(), yNeg->GetPy(), yNeg->GetPz());
2408
2409     Double_t deltaTheta = neg.Theta() - pos.Theta();
2410     Double_t openingAngle = pos.Angle( neg );
2411
2412     if( openingAngle < 1e-20 ) return 0.;
2413
2414     return TMath::ASin( deltaTheta/openingAngle );
2415 }
2416
2417 //-----------------------------------------------------------------------------------------------
2418 Double_t AliAnalysisTaskGammaConvDalitz::GetPsiPair(const TLorentzVector* xPos, const TLorentzVector* yNeg ) const
2419 {
2420 //
2421 // Get psi pair value
2422 //
2423         Double_t deltaTheta = yNeg->Theta() - xPos->Theta();
2424         Double_t openingAngle = xPos->Angle( yNeg->Vect() );
2425
2426         if( openingAngle < 1e-20 ) return 0.;
2427
2428         return TMath::ASin( deltaTheta/openingAngle );;
2429 }
2430
2431 // tmp  NOTE: Should go to AliV0Reader
2432 //-----------------------------------------------------------------------------------------------
2433 /*
2434 Double_t AliAnalysisTaskGammaConvDalitz::GetPsiAngleV0(
2435             Double_t  radiusVO,                     // radius at XY of VO vertex
2436             const AliExternalTrackParam* trackPos,  // pos. track parm. at V0 vertex
2437             const AliExternalTrackParam* trackNeg  // neg. track parm. at V0 vertex
2438          )
2439 {
2440     // This angle is a measure for the contribution of the opening in polar
2441     // direction Δ0 to the opening angle ξ Pair
2442
2443     // Ref. Measurement of photons via conversion pairs with the PHENIX experiment at RHIC
2444     //      Master Thesis. Thorsten Dahms. 2005
2445     // https://twiki.cern.ch/twiki/pub/ALICE/GammaPhysicsPublications/tdahms_thesis.pdf
2446
2447     const Double_t kForceDeltaPhi = 1.2;
2448     static const Double_t kB2C = 0.299792458e-3; // Taken from AliVParticle
2449
2450     Double_t psiAngle = 0.;  // Default value
2451
2452     static Double_t MagnField = fESDEvent->GetMagneticField();
2453     static Double_t MagnFieldG = fESDEvent->GetMagneticField()*kB2C;
2454
2455     if( TMath::Abs( MagnField ) < 1e-20 ) return psiAngle; // Do nothing if 0 mag field
2456
2457     // compute propagation radius for a fixed angle
2458     Double_t Rpos = radiusVO + trackPos->Pt()*TMath::Sin( kForceDeltaPhi ) / MagnFieldG;
2459     Double_t Rneg = radiusVO + trackNeg->Pt()*TMath::Sin( kForceDeltaPhi ) / MagnFieldG;
2460
2461     Double_t MomPos[3];
2462     Double_t MomNeg[3];
2463     if( trackPos->GetPxPyPzAt( Rpos, MagnField, MomPos ) == 0 ) trackPos->GetPxPyPz( MomPos );
2464     if( trackNeg->GetPxPyPzAt( Rneg, MagnField, MomNeg ) == 0 ) trackNeg->GetPxPyPz( MomNeg );
2465
2466     TVector3 PosDaughter;
2467     TVector3 NegDaughter;
2468     PosDaughter.SetXYZ( MomPos[0], MomPos[1], MomPos[2] );
2469     NegDaughter.SetXYZ( MomNeg[0], MomNeg[1], MomNeg[2] );
2470
2471     Double_t deltaTheta = NegDaughter.Theta() - PosDaughter.Theta();
2472     Double_t chiPar = PosDaughter.Angle( NegDaughter );  //TMath::ACos( PosDaughter.Dot(NegDaughter) / (NegDaughter.Mag()*PosDaughter.Mag()) );
2473     if( chiPar < 1e-20 ) return psiAngle;
2474
2475     psiAngle = TMath::ASin( deltaTheta / chiPar );
2476
2477     return psiAngle;
2478 }
2479 */