PWGDQ/dielectron/AliDielectronMC.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 //#####################################################
  17 //#                                                   #
  18 //#              Class AliDielectronMC                #
  19 //#       Cut Class for Jpsi->e+e- analysis           #
  20 //#                                                   #
  21 //#   by WooJin J. Park, GSI / W.J.Park@gsi.de        #
  22 //#                                                   #
  23 //#####################################################
  24
  25 #include <AliAnalysisManager.h>
  26 #include <AliAODHandler.h>
  27 #include <AliESDInputHandler.h>
  28 #include <AliAODInputHandler.h>
  29 #include <AliMCEventHandler.h>
  30 #include <AliMCEvent.h>
  31 #include <AliMCParticle.h>
  32 #include <AliAODMCParticle.h>
  33 #include <AliAODMCHeader.h>
  34 #include <AliStack.h>
  35 #include <AliESDEvent.h>
  36 #include <AliAODEvent.h>
  37 #include <AliESDtrack.h>
  38 #include <AliAODTrack.h>
  39 #include <AliLog.h>
  40
  41 #include <TClonesArray.h>
  42 #include <TParticle.h>
  43
  44 #include "AliDielectronSignalMC.h"
  45 #include "AliDielectronMC.h"
  46
  47
  48 ClassImp(AliDielectronMC)
  49
  50 AliDielectronMC* AliDielectronMC::fgInstance=0x0;
  51
  52 //____________________________________________________________
  53 AliDielectronMC* AliDielectronMC::Instance()
  54 {
  55   //
  56   // return pointer to singleton implementation
  57   //
  58   if (fgInstance) return fgInstance;
  59
  60   AnalysisType type=kUNSET;
  61   Bool_t hasMC=kFALSE;
  62   if (AliAnalysisManager::GetAnalysisManager()){
  63     if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliESDInputHandler::Class()) type=kESD;
  64     else if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliAODInputHandler::Class()) type=kAOD;
  65
  66     AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
  67     if(type == kESD) hasMC=mcHandler!=0x0;
  68     }
  69
  70   fgInstance=new AliDielectronMC(type);
  71
  72   fgInstance->SetHasMC(hasMC);
  73
  74   return fgInstance;
  75 }
  76
  77 //____________________________________________________________
  78 AliDielectronMC::AliDielectronMC(AnalysisType type):
  79   fMCEvent(0x0),
  80   fStack(0x0),
  81   fAnaType(type),
  82   fHasMC(kTRUE),
  83   fMcArray(0x0)
  84 {
  85   //
  86   // default constructor
  87   //
  88 }
  89
  90
  91 //____________________________________________________________
  92 AliDielectronMC::~AliDielectronMC()
  93 {
  94   //
  95   // default destructor
  96   //
  97
  98 }
  99
 100 //____________________________________________________________
 101 void AliDielectronMC::Initialize()
 102 {
 103   //
 104   // initialize MC class
 105   //
 106   if (!ConnectMCEvent()) AliError("Initialization of MC object failed!");
 107 }
 108
 109 //____________________________________________________________
 110 Int_t AliDielectronMC::GetNMCTracks()
 111 {
 112   //
 113   //  return the number of generated tracks from MC event
 114   //
 115   if(fAnaType == kESD){
 116     if (!fMCEvent){ AliError("No fMCEvent"); return 0; }
 117     return fMCEvent->GetNumberOfTracks();}
 118   else if(fAnaType == kAOD){
 119     if(!fMcArray) { AliError("No fMcArray"); return 0; }
 120     return fMcArray->GetEntriesFast();
 121   }
 122   return 0;
 123 }
 124
 125 //____________________________________________________________
 126 Int_t AliDielectronMC::GetNMCTracksFromStack()
 127 {
 128   //
 129   //  return the number of generated tracks from stack
 130   //
 131   if (!fStack){ AliError("No fStack"); return -999; }
 132   return fStack->GetNtrack();
 133 }
 134
 135 //____________________________________________________________
 136 Int_t AliDielectronMC::GetNPrimary() const
 137 {
 138   //
 139   //  return the number of primary track from MC event
 140   //
 141   if (!fMCEvent){ AliError("No fMCEvent"); return 0; }
 142   return fMCEvent->GetNumberOfPrimaries();
 143 }
 144
 145 //____________________________________________________________
 146 Int_t AliDielectronMC::GetNPrimaryFromStack()
 147 {
 148   //
 149   //  return the number of primary track from stack
 150   //
 151   if (!fStack){ AliError("No fStack"); return -999; }
 152   return fStack->GetNprimary();
 153 }
 154
 155 //____________________________________________________________
 156 AliVParticle* AliDielectronMC::GetMCTrackFromMCEvent(Int_t label) const
 157 {
 158   //
 159   // return MC track directly from MC event
 160   // used not only for tracks but for mothers as well, therefore do not use abs(label)
 161   //
 162   if (label<0) return NULL;
 163   AliVParticle * track=0x0;
 164   if(fAnaType == kESD){
 165     if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
 166     track = fMCEvent->GetTrack(label); //  tracks from MC event (ESD)
 167   } else if(fAnaType == kAOD) {
 168     if (!fMcArray){ AliError("No fMcArray"); return NULL;}
 169     if (label>fMcArray->GetEntriesFast()) { AliDebug(10,Form("track %d out of array size %d",label,fMcArray->GetEntriesFast())); return NULL;}
 170     track = (AliVParticle*)fMcArray->At(label); //  tracks from MC event (AOD)
 171   }
 172   return track;
 173 }
 174
 175 //____________________________________________________________
 176 Bool_t AliDielectronMC::ConnectMCEvent()
 177 {
 178   //
 179   // connect stack object from the mc handler
 180   //
 181
 182   fMcArray = 0x0;
 183   fMCEvent = 0x0;
 184
 185   if(fAnaType == kESD){
 186     AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 187     if (!mcHandler){ /*AliError("Could not retrive MC event handler!");*/ return kFALSE; }
 188     if (!mcHandler->InitOk() ) return kFALSE;
 189     if (!mcHandler->TreeK() )  return kFALSE;
 190     if (!mcHandler->TreeTR() ) return kFALSE;
 191
 192     AliMCEvent* mcEvent = mcHandler->MCEvent();
 193     if (!mcEvent){ /*AliError("Could not retrieve MC event!");*/ return kFALSE; }
 194     fMCEvent = mcEvent;
 195
 196     if (!UpdateStack()) return kFALSE;
 197   }
 198   else if(fAnaType == kAOD)
 199   {
 200     AliAODInputHandler* aodHandler=(AliAODInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 201     if (!aodHandler) return kFALSE;
 202     AliAODEvent *aod=aodHandler->GetEvent();
 203     if (!aod) return kFALSE;
 204
 205     fMcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
 206     if (!fMcArray){ /*AliError("Could not retrieve MC array!");*/ return kFALSE; }
 207     else fHasMC=kTRUE;
 208   }
 209   return kTRUE;
 210 }
 211
 212 //____________________________________________________________
 213 Bool_t AliDielectronMC::UpdateStack()
 214 {
 215   //
 216   // update stack with new event
 217   //
 218   if (!fMCEvent){ AliError("No fMCEvent"); return kFALSE;}
 219   AliStack* stack = fMCEvent->Stack();
 220   if (!stack){ AliError("Could not retrive stack!"); return kFALSE; }
 221   fStack = stack;
 222   return kTRUE;
 223 }
 224
 225 //____________________________________________________________
 226 AliMCParticle* AliDielectronMC::GetMCTrack( const AliESDtrack* _track)
 227 {
 228   //
 229   // return MC track
 230   //
 231   if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
 232   Int_t nStack = fMCEvent->GetNumberOfTracks();
 233   Int_t label  = TMath::Abs(_track->GetLabel()); // negative label indicate poor matching quality
 234   if(label>nStack)return NULL;
 235   AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(label));
 236   return mctrack;
 237 }
 238
 239
 240 //____________________________________________________________
 241 AliAODMCParticle* AliDielectronMC::GetMCTrack( const AliAODTrack* _track)
 242 {
 243   //
 244   // return MC track
 245   //
 246  if(!fMcArray) { AliError("No fMcArray"); return NULL;}
 247  Int_t nStack = fMcArray->GetEntriesFast();
 248  Int_t label  = TMath::Abs(_track->GetLabel()); // negative label indicate poor matching quality
 249  if(label > nStack) return NULL;
 250  AliAODMCParticle *mctrack = (AliAODMCParticle*)fMcArray->At(label);
 251  return mctrack;
 252 }
 253
 254 //____________________________________________________________
 255 TParticle* AliDielectronMC::GetMCTrackFromStack(const AliESDtrack* _track)
 256 {
 257   //
 258   // return MC track from stack
 259   //
 260   Int_t label = TMath::Abs(_track->GetLabel());
 261   if (!fStack) AliWarning("fStack is not available. Update stack first.");
 262   TParticle* mcpart = fStack->Particle(label);
 263   if (!mcpart) return NULL;
 264   return mcpart;
 265 }
 266
 267 //____________________________________________________________
 268 AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliESDtrack* _track)
 269 {
 270   //
 271   // return MC track mother
 272   //
 273   AliMCParticle* mcpart = GetMCTrack(_track);
 274   if (!mcpart) return NULL;
 275   if(mcpart->GetMother()<0) return NULL;
 276   AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(mcpart->GetMother()));
 277   if (!mcmother) return NULL;
 278   return mcmother;
 279 }
 280
 281 //______________________________________________________________
 282 AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODTrack* _track)
 283 {
 284  //
 285  // return MC track mother
 286  //
 287  AliAODMCParticle* mcpart = GetMCTrack(_track);
 288  if (!mcpart) return NULL;
 289  if(mcpart->GetMother() < 0) return NULL;
 290  AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMcArray->At(mcpart->GetMother()));
 291  if (!mcmother) return NULL;
 292  return mcmother;
 293 }
 294 //____________________________________________________________
 295 AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliMCParticle* _particle){
 296   //
 297   // return MC track mother
 298   //
 299   if(_particle->GetMother() < 0) return NULL;
 300   AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(_particle->GetMother()));
 301   return mcmother;
 302 }
 303
 304 //____________________________________________________________
 305 AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODMCParticle* _particle){
 306   //
 307   // return MC track mother
 308   //
 309   if( _particle->GetMother() < 0) return NULL;
 310   AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMcArray->At(_particle->GetMother()));
 311   return mcmother;
 312 }
 313
 314 //____________________________________________________________
 315 TParticle* AliDielectronMC::GetMCTrackMotherFromStack(const AliESDtrack* _track)
 316 {
 317   //
 318   // return MC track mother from stack
 319   //
 320   TParticle* mcpart = GetMCTrackFromStack(_track);
 321   if ( !mcpart || mcpart->GetFirstMother()<=0 ) return NULL;
 322   TParticle* mcmother = fStack->Particle(mcpart->GetFirstMother());
 323   if (!mcmother) return NULL;
 324   return mcmother;
 325 }
 326
 327 //____________________________________________________________
 328 Int_t AliDielectronMC::GetMCPID(const AliESDtrack* _track)
 329 {
 330   //
 331   // return PDG code of the track from the MC truth info
 332   //
 333   AliMCParticle* mcpart = GetMCTrack(_track);
 334   if (!mcpart) return -999;
 335   return mcpart->PdgCode();
 336 }
 337
 338 //__________________________________________________________
 339 Int_t AliDielectronMC::GetMCPID(const AliAODTrack* _track)
 340 {
 341  //
 342  // return PDG code of the track from the MC truth info
 343  //
 344  AliAODMCParticle* mcpart = GetMCTrack(_track);
 345  if (!mcpart) return -999;
 346  return mcpart->PdgCode();
 347 }
 348
 349 //____________________________________________________________
 350 Int_t AliDielectronMC::GetMCPIDFromStack(const AliESDtrack* _track)
 351 {
 352   //
 353   // return MC PDG code from stack
 354   //
 355   TParticle* mcpart = GetMCTrackFromStack(_track);
 356   if (!mcpart) return -999;
 357   return mcpart->GetPdgCode();
 358 }
 359
 360 //____________________________________________________________
 361 Int_t AliDielectronMC::GetMotherPDG( const AliESDtrack* _track)
 362 {
 363   //
 364   // return PDG code of the mother track from the MC truth info
 365   //
 366   AliMCParticle* mcmother = GetMCTrackMother(_track);
 367   if (!mcmother) return -999;
 368   return mcmother->PdgCode();
 369 }
 370
 371 //________________________________________________________
 372 Int_t AliDielectronMC::GetMotherPDG( const AliAODTrack* _track)
 373 {
 374   //
 375   // return PDG code of the mother track from the MC truth info
 376   //
 377   AliAODMCParticle* mcmother = GetMCTrackMother(_track);
 378   if (!mcmother) return -999;
 379   return mcmother->PdgCode();
 380 }
 381
 382 //________________________________________________________
 383 Int_t AliDielectronMC::GetMotherPDG( const AliMCParticle* _track)
 384 {
 385   //
 386   // return PDG code of the mother track from the MC truth info
 387   //
 388   AliMCParticle* mcmother = GetMCTrackMother(_track);
 389   if (!mcmother) return -999;
 390   return mcmother->PdgCode();
 391 }
 392
 393 //________________________________________________________
 394 Int_t AliDielectronMC::GetMotherPDG( const AliAODMCParticle* _track)
 395 {
 396   //
 397   // return PDG code of the mother track from the MC truth info
 398   //
 399   AliAODMCParticle* mcmother = GetMCTrackMother(_track);
 400   if (!mcmother) return -999;
 401   return mcmother->PdgCode();
 402 }
 403
 404 //____________________________________________________________
 405 Int_t AliDielectronMC::GetMotherPDGFromStack(const AliESDtrack* _track)
 406 {
 407   //
 408   // return PDG code of the mother track from stack
 409   //
 410   TParticle* mcmother = GetMCTrackMotherFromStack(_track);
 411   if (!mcmother) return -999;
 412   return mcmother->GetPdgCode();
 413 }
 414
 415 //____________________________________________________________
 416 Int_t AliDielectronMC::GetMCProcess(const AliESDtrack* _track)
 417 {
 418   //
 419   // return process number of the track
 420   //
 421   AliMCParticle* mcpart = GetMCTrack(_track);
 422   if (!mcpart) return -999;
 423   return 0;
 424 }
 425
 426 //____________________________________________________________
 427 Int_t AliDielectronMC::GetMCProcessFromStack(const AliESDtrack* _track)
 428 {
 429   //
 430   // return process number of the track
 431   //
 432   TParticle* mcpart = GetMCTrackFromStack(_track);
 433   if (!mcpart) return -999;
 434   return mcpart->GetUniqueID();
 435 }
 436
 437 //____________________________________________________________
 438 Int_t AliDielectronMC::NumberOfDaughters(const AliESDtrack* track)
 439 {
 440   //
 441   // returns the number of daughters
 442   //
 443   AliMCParticle *mcmother=GetMCTrackMother(track);
 444   if(!mcmother||!mcmother->Particle()) return -999;
 445   //   return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
 446   return mcmother->Particle()->GetNDaughters();
 447 }
 448
 449 //_________________________________________________________
 450 Int_t AliDielectronMC::NumberOfDaughters(const AliAODTrack* track)
 451 {
 452   //
 453   // returns the number of daughters
 454   //
 455   AliAODMCParticle *mcmother=GetMCTrackMother(track);
 456   if(!mcmother) return -999;
 457   return NumberOfDaughters(mcmother);
 458
 459 }
 460
 461 //____________________________________________________________
 462 Int_t AliDielectronMC::NumberOfDaughters(const AliMCParticle* particle)
 463 {
 464   //
 465   // returns the number of daughters
 466   //
 467   AliMCParticle *mcmother=GetMCTrackMother(particle);
 468   if(!mcmother||!mcmother->Particle()) return -999;
 469   //return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
 470   return mcmother->Particle()->GetNDaughters();
 471 }
 472
 473 //____________________________________________________________
 474 Int_t AliDielectronMC::NumberOfDaughters(const AliAODMCParticle* particle)
 475 {
 476   //
 477   // returns the number of daughters
 478   //
 479   AliAODMCParticle *mcmother=GetMCTrackMother(particle);
 480   if(!mcmother) return -999;
 481   return mcmother->GetNDaughters();
 482 }
 483
 484 //____________________________________________________________
 485 Int_t AliDielectronMC::GetMCProcessMother(const AliESDtrack* _track)
 486 {
 487   //
 488   // return process number of the mother of the track
 489   //
 490   AliMCParticle* mcmother = GetMCTrackMother(_track);
 491   if (!mcmother) return -999;
 492   return 0;
 493 }
 494
 495 //____________________________________________________________
 496 Int_t AliDielectronMC::GetMCProcessMotherFromStack(const AliESDtrack* _track)
 497 {
 498   //
 499   // return process number of the mother of the track
 500   //
 501   TParticle* mcmother = GetMCTrackMotherFromStack(_track);
 502   if (!mcmother) return -999;
 503   return mcmother->GetUniqueID();
 504 }
 505
 506 //____________________________________________________________
 507 Bool_t AliDielectronMC::IsMCMotherToEE(const AliVParticle *particle, Int_t pdgMother)
 508 {
 509   //
 510   // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
 511   //
 512   if (fAnaType==kESD && !fMCEvent) return kFALSE;
 513   if (fAnaType==kAOD && !fMcArray) return kFALSE;
 514   if (!particle) return kFALSE;
 515
 516   if (particle->IsA()==AliMCParticle::Class()){
 517     return IsMCMotherToEEesd(static_cast<const AliMCParticle*>(particle),pdgMother);
 518   } else if (particle->IsA()==AliAODMCParticle::Class()){
 519    return IsMCMotherToEEaod(static_cast<const AliAODMCParticle*>(particle),pdgMother);
 520   } else {
 521     AliError("Unknown particle type");
 522   }
 523   return kFALSE;
 524 }
 525
 526 //____________________________________________________________
 527 Bool_t AliDielectronMC::IsMCMotherToEEesd(const AliMCParticle *particle, Int_t pdgMother)
 528 {
 529   //
 530   // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
 531   // ESD case
 532   //
 533
 534   //check pdg code
 535   if (particle->PdgCode()!=pdgMother) return kFALSE;
 536   Int_t ifirst = particle->GetFirstDaughter();
 537   Int_t ilast  = particle->GetLastDaughter();
 538
 539   //check number of daughters
 540   if ((ilast-ifirst)!=1) return kFALSE;
 541   AliMCParticle *firstD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ifirst));
 542   AliMCParticle *secondD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ilast));
 543
 544   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
 545   if (firstD->Charge()>0){
 546     if (firstD->PdgCode()!=-11) return kFALSE;
 547     if (secondD->PdgCode()!=11) return kFALSE;
 548   }else{
 549     if (firstD->PdgCode()!=11) return kFALSE;
 550     if (secondD->PdgCode()!=-11) return kFALSE;
 551   }
 552
 553   return kTRUE;
 554 }
 555
 556 //____________________________________________________________
 557 Bool_t AliDielectronMC::IsMCMotherToEEaod(const AliAODMCParticle *particle, Int_t pdgMother)
 558 {
 559   //
 560   // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
 561   // AOD case
 562   //
 563
 564   if (particle->GetPdgCode()!=pdgMother) return kFALSE;
 565   if (particle->GetNDaughters()!=2) return kFALSE;
 566
 567   Int_t ifirst = particle->GetDaughter(0);
 568   Int_t ilast  = particle->GetDaughter(1);
 569
 570   //check number of daughters
 571   if ((ilast-ifirst)!=1) return kFALSE;
 572
 573   AliAODMCParticle *firstD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ifirst));
 574   AliAODMCParticle *secondD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ilast));
 575
 576   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
 577
 578   if (firstD->Charge()>0){
 579     if (firstD->GetPdgCode()!=-11) return kFALSE;
 580     if (secondD->GetPdgCode()!=11) return kFALSE;
 581   }else{
 582     if (firstD->GetPdgCode()!=11) return kFALSE;
 583     if (secondD->GetPdgCode()!=-11) return kFALSE;
 584   }
 585   return kTRUE;
 586 }
 587
 588 //____________________________________________________________
 589 Int_t AliDielectronMC::GetLabelMotherWithPdg(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
 590 {
 591   //
 592   // test if mother of particle 1 and 2 has pdgCode pdgMother and is the same;
 593   //
 594   if (fAnaType==kESD){
 595   if (!fMCEvent) return -1;
 596   return GetLabelMotherWithPdgESD(particle1, particle2, pdgMother);
 597   }
 598   else if (fAnaType==kAOD)
 599   {
 600   if (!fMcArray) return -1;
 601   return GetLabelMotherWithPdgAOD(particle1, particle2, pdgMother);
 602   }
 603
 604   return -1;
 605 }
 606
 607 //____________________________________________________________
 608 Int_t AliDielectronMC::GetLabelMotherWithPdgESD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
 609 {
 610   //
 611   // test if mother of particle 1 and 2 has pdgCode +-11 (electron),
 612   //    have the same mother and the mother had pdg code pdgMother
 613   // ESD case
 614   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
 615   //
 616   // negative label indicate poor matching quality
 617   Int_t lblPart1 = TMath::Abs(particle1->GetLabel());
 618   Int_t lblPart2 = TMath::Abs(particle2->GetLabel());
 619   AliMCParticle *mcPart1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(lblPart1));
 620   AliMCParticle *mcPart2=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(lblPart2));
 621
 622   if (!mcPart1||!mcPart2) return -1;
 623
 624   Int_t lblMother1=mcPart1->GetMother();
 625   Int_t lblMother2=mcPart2->GetMother();
 626   AliMCParticle *mcMother1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(lblMother1));
 627
 628   if (!mcMother1) return -1;
 629   if (lblMother1!=lblMother2) return -1;
 630   if (TMath::Abs(mcPart1->PdgCode())!=11) return -1;
 631   if (mcPart1->PdgCode()!=-mcPart2->PdgCode()) return -1;
 632   if (mcMother1->PdgCode()!=pdgMother) return -1;
 633
 634   return lblMother1;
 635 }
 636
 637 //____________________________________________________________
 638 Int_t AliDielectronMC::GetLabelMotherWithPdgAOD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
 639 {
 640   //
 641   // test if mother of particle 1 and 2 has pdgCode +-11 (electron),
 642   //    have the same mother and the mother had pdg code pdgMother
 643   // AOD case
 644   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
 645   //
 646   // negative label indicate poor matching quality
 647   Int_t lblPart1 = TMath::Abs(particle1->GetLabel());
 648   Int_t lblPart2 = TMath::Abs(particle2->GetLabel());
 649   AliAODMCParticle *mcPart1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(lblPart1));
 650   AliAODMCParticle *mcPart2=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(lblPart2));
 651
 652   if (!mcPart1||!mcPart2) return -1;
 653
 654   Int_t lblMother1=mcPart1->GetMother();
 655   Int_t lblMother2=mcPart2->GetMother();
 656   AliAODMCParticle *mcMother1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(lblMother1));
 657
 658   if (!mcMother1) return -1;
 659   if (lblMother1!=lblMother2) return -1;
 660   if (TMath::Abs(mcPart1->GetPdgCode())!=11) return -1;
 661   if (mcPart1->GetPdgCode()!=-mcPart2->GetPdgCode()) return -1;
 662   if (mcMother1->GetPdgCode()!=pdgMother) return -1;
 663
 664   return lblMother1;
 665 }
 666
 667 //____________________________________________________________
 668 void AliDielectronMC::GetDaughters(const TObject *mother, AliVParticle* &d1, AliVParticle* &d2)
 669 {
 670   //
 671   // Get First two daughters of the mother
 672   //
 673   Int_t lblD1=-1;
 674   Int_t lblD2=-1;
 675   d1=0;
 676   d2=0;
 677   if (fAnaType==kAOD){
 678     if(!fMcArray) return;
 679     const AliAODMCParticle *aodMother=static_cast<const AliAODMCParticle*>(mother);
 680     lblD1=aodMother->GetDaughter(0);
 681     lblD2=aodMother->GetDaughter(1);
 682     d1 = (AliVParticle*)fMcArray->At(lblD1);
 683     d2 = (AliVParticle*)fMcArray->At(lblD2);
 684    } else if (fAnaType==kESD){
 685     if (!fMCEvent) return;
 686     const AliMCParticle *aodMother=static_cast<const AliMCParticle*>(mother);
 687     lblD1=aodMother->GetFirstDaughter();
 688     lblD2=aodMother->GetLastDaughter();
 689     d1=fMCEvent->GetTrack(lblD1);
 690     d2=fMCEvent->GetTrack(lblD2);
 691    }
 692 }
 693
 694
 695 //________________________________________________________________________________
 696 Int_t AliDielectronMC::GetMothersLabel(Int_t daughterLabel) const {
 697   //
 698   //  Get the label of the mother for particle with label daughterLabel
 699   //  NOTE: for tracks, the absolute label should be passed
 700   //
 701   if(daughterLabel<0) return -1;
 702   if (fAnaType==kAOD) {
 703     if(!fMcArray) return -1;
 704     if(GetMCTrackFromMCEvent(daughterLabel))
 705       return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(daughterLabel)))->GetMother();
 706   } else if(fAnaType==kESD) {
 707     if (!fMCEvent) return -1;
 708     if(GetMCTrackFromMCEvent(daughterLabel))
 709       return (static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(daughterLabel)))->GetMother();
 710   }
 711   return -1;
 712 }
 713
 714
 715 //________________________________________________________________________________
 716 Int_t AliDielectronMC::GetPdgFromLabel(Int_t label) const {
 717   //
 718   //  Get particle code using the label from stack
 719   //  NOTE: for tracks, the absolute label should be passed
 720   //
 721   if(label<0) return 0;
 722   if(fAnaType==kAOD) {
 723     if(!fMcArray) return 0;
 724     return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label)))->PdgCode();
 725   } else if(fAnaType==kESD) {
 726     if (!fMCEvent) return 0;
 727     return (static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(label)))->PdgCode();
 728   }
 729   return 0;
 730 }
 731
 732
 733 //________________________________________________________________________________
 734 Bool_t AliDielectronMC::ComparePDG(Int_t particlePDG, Int_t requiredPDG, Bool_t pdgExclusion, Bool_t checkBothCharges) const {
 735   //
 736   //  Test the PDG codes of particles with the required ones
 737   //
 738   Bool_t result = kTRUE;
 739   Int_t absRequiredPDG = TMath::Abs(requiredPDG);
 740
 741   switch(absRequiredPDG) {
 742   case 0:
 743     result = kTRUE;    // PDG not required (any code will do fine)
 744     break;
 745   case 100:     // light flavoured mesons
 746     if(checkBothCharges)
 747       result = TMath::Abs(particlePDG)>=100 && TMath::Abs(particlePDG)<=199;
 748     else {
 749       if(requiredPDG>0) result = particlePDG>=100 && particlePDG<=199;
 750       if(requiredPDG<0) result = particlePDG>=-199 && particlePDG<=-100;
 751     }
 752     break;
 753   case 1000:     // light flavoured baryons
 754     if(checkBothCharges)
 755       result = TMath::Abs(particlePDG)>=1000 && TMath::Abs(particlePDG)<=1999;
 756     else {
 757       if(requiredPDG>0) result = particlePDG>=1000 && particlePDG<=1999;
 758       if(requiredPDG<0) result = particlePDG>=-1999 && particlePDG<=-1000;
 759     }
 760     break;
 761   case 200:     // light flavoured mesons
 762     if(checkBothCharges)
 763       result = TMath::Abs(particlePDG)>=200 && TMath::Abs(particlePDG)<=299;
 764     else {
 765       if(requiredPDG>0)result = particlePDG>=200 && particlePDG<=299;
 766       if(requiredPDG<0)result = particlePDG>=-299 && particlePDG<=-200;
 767     }
 768     break;
 769   case 2000:     // light flavoured baryons
 770     if(checkBothCharges)
 771       result = TMath::Abs(particlePDG)>=2000 && TMath::Abs(particlePDG)<=2999;
 772     else {
 773       if(requiredPDG>0) result = particlePDG>=2000 && particlePDG<=2999;
 774       if(requiredPDG<0) result = particlePDG>=-2999 && particlePDG<=-2000;
 775     }
 776     break;
 777   case 300:     // all strange mesons
 778     if(checkBothCharges)
 779       result = TMath::Abs(particlePDG)>=300 && TMath::Abs(particlePDG)<=399;
 780     else {
 781       if(requiredPDG>0) result = particlePDG>=300 && particlePDG<=399;
 782       if(requiredPDG<0) result = particlePDG>=-399 && particlePDG<=-300;
 783     }
 784     break;
 785   case 3000:     // all strange baryons
 786     if(checkBothCharges)
 787       result = TMath::Abs(particlePDG)>=3000 && TMath::Abs(particlePDG)<=3999;
 788     else {
 789       if(requiredPDG>0) result = particlePDG>=3000 && particlePDG<=3999;
 790       if(requiredPDG<0) result = particlePDG>=-3999 && particlePDG<=-3000;
 791     }
 792     break;
 793   case 400:     // all charmed mesons
 794     if(checkBothCharges)
 795       result = TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=499;
 796     else {
 797       if(requiredPDG>0) result = particlePDG>=400 && particlePDG<=499;
 798       if(requiredPDG<0) result = particlePDG>=-499 && particlePDG<=-400;
 799     }
 800     break;
 801   case 401:     // open charm mesons
 802     if(checkBothCharges)
 803       result = TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=439;
 804     else {
 805       if(requiredPDG>0) result = particlePDG>=400 && particlePDG<=439;
 806       if(requiredPDG<0) result = particlePDG>=-439 && particlePDG<=-400;
 807     }
 808     break;
 809   case 402:     // open charm mesons and baryons together
 810     if(checkBothCharges)
 811       result = (TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=439) ||
 812                (TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4399);
 813     else {
 814       if(requiredPDG>0) result = (particlePDG>=400 && particlePDG<=439) ||
 815                                  (particlePDG>=4000 && particlePDG<=4399);
 816       if(requiredPDG<0) result = (particlePDG>=-439 && particlePDG<=-400) ||
 817                                  (particlePDG>=-4399 && particlePDG<=-4000);
 818     }
 819     break;
 820   case 403:     // all charm hadrons
 821     if(checkBothCharges)
 822       result = (TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=499) ||
 823                (TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4999);
 824     else {
 825       if(requiredPDG>0) result = (particlePDG>=400 && particlePDG<=499) ||
 826                                  (particlePDG>=4000 && particlePDG<=4999);
 827       if(requiredPDG<0) result = (particlePDG>=-499 && particlePDG<=-400) ||
 828                                  (particlePDG>=-4999 && particlePDG<=-4000);
 829     }
 830     break;
 831   case 4000:     // all charmed baryons
 832     if(checkBothCharges)
 833       result = TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4999;
 834     else {
 835       if(requiredPDG>0) result = particlePDG>=4000 && particlePDG<=4999;
 836       if(requiredPDG<0) result = particlePDG>=-4999 && particlePDG<=-4000;
 837     }
 838     break;
 839   case 4001:     // open charm baryons
 840     if(checkBothCharges)
 841       result = TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4399;
 842     else {
 843       if(requiredPDG>0) result = particlePDG>=4000 && particlePDG<=4399;
 844       if(requiredPDG<0) result = particlePDG>=-4399 && particlePDG<=-4000;
 845     }
 846     break;
 847   case 500:      // all beauty mesons
 848     if(checkBothCharges)
 849       result = TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=599;
 850     else {
 851       if(requiredPDG>0) result = particlePDG>=500 && particlePDG<=599;
 852       if(requiredPDG<0) result = particlePDG>=-599 && particlePDG<=-500;
 853     }
 854     break;
 855   case 501:      // open beauty mesons
 856     if(checkBothCharges)
 857       result = TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=549;
 858     else {
 859       if(requiredPDG>0) result = particlePDG>=500 && particlePDG<=549;
 860       if(requiredPDG<0) result = particlePDG>=-549 && particlePDG<=-500;
 861     }
 862     break;
 863   case 502:      // open beauty mesons and baryons
 864     if(checkBothCharges)
 865       result = (TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=549) ||
 866                (TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5499);
 867     else {
 868       if(requiredPDG>0) result = (particlePDG>=500 && particlePDG<=549) ||
 869                                  (particlePDG>=5000 && particlePDG<=5499);
 870       if(requiredPDG<0) result = (particlePDG>=-549 && particlePDG<=-500) ||
 871                                  (particlePDG>=-5499 && particlePDG<=-5000);
 872     }
 873     break;
 874   case 503:      // all beauty hadrons
 875     if(checkBothCharges)
 876       result = (TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=599) ||
 877                (TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5999);
 878     else {
 879       if(requiredPDG>0) result = (particlePDG>=500 && particlePDG<=599) ||
 880                                  (particlePDG>=5000 && particlePDG<=5999);
 881       if(requiredPDG<0) result = (particlePDG>=-599 && particlePDG<=-500) ||
 882                                  (particlePDG>=-5999 && particlePDG<=-5000);
 883     }
 884     break;
 885   case 5000:      // all beauty baryons
 886     if(checkBothCharges)
 887       result = TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5999;
 888     else {
 889       if(requiredPDG>0) result = particlePDG>=5000 && particlePDG<=5999;
 890       if(requiredPDG<0) result = particlePDG>=-5999 && particlePDG<=-5000;
 891     }
 892     break;
 893   case 5001:      // open beauty baryons
 894     if(checkBothCharges)
 895       result = TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5499;
 896     else {
 897       if(requiredPDG>0) result = particlePDG>=5000 && particlePDG<=5499;
 898       if(requiredPDG<0) result = particlePDG>=-5499 && particlePDG<=-5000;
 899     }
 900     break;
 901   case 902:      // // open charm,beauty  mesons and baryons together
 902     if(checkBothCharges)
 903       result = (TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=439) ||
 904         (TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4399) ||
 905         (TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=549) ||
 906         (TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5499);
 907     else {
 908       if(requiredPDG>0) result = (particlePDG>=400 && particlePDG<=439) ||
 909                           (particlePDG>=4000 && particlePDG<=4399)      ||
 910                           (particlePDG>=500 && particlePDG<=549)        ||
 911                           (particlePDG>=5000 && particlePDG<=5499);
 912       if(requiredPDG<0) result = (particlePDG>=-439 && particlePDG<=-400) ||
 913                           (particlePDG>=-4399 && particlePDG<=-4000)      ||
 914                           (particlePDG>=-549 && particlePDG<=-500)        ||
 915                           (particlePDG>=-5499 && particlePDG<=-5000);
 916     }
 917     break;
 918   default:          // all specific cases
 919     if(checkBothCharges)
 920       result = (absRequiredPDG==TMath::Abs(particlePDG));
 921     else
 922       result = (requiredPDG==particlePDG);
 923   }
 924
 925   if(absRequiredPDG!=0 && pdgExclusion) result = !result;
 926   return result;
 927 }
 928
 929 //________________________________________________________________________________
 930 Bool_t AliDielectronMC::IsPhysicalPrimary(Int_t label) const {
 931   //
 932   // Check if the particle with label "label" is a physical primary according to the
 933   // definition in AliStack::IsPhysicalPrimary(Int_t label)
 934   // Convention for being physical primary:
 935   // 1.) particles produced in the collision
 936   // 2.) stable particles with respect to strong and electromagnetic interactions
 937   // 3.) excludes initial state particles
 938   // 4.) includes products of directly produced Sigma0 hyperon decay
 939   // 5.) includes products of directly produced pi0 decays
 940   // 6.) includes products of directly produced beauty hadron decays
 941   //
 942   if(label<0) return kFALSE;
 943   if(fAnaType==kAOD) {
 944     if(!fMcArray) return kFALSE;
 945     return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label)))->IsPhysicalPrimary();
 946   } else if(fAnaType==kESD) {
 947     if (!fMCEvent) return kFALSE;
 948     return fStack->IsPhysicalPrimary(label);
 949   }
 950   return kFALSE;
 951 }
 952
 953
 954 //________________________________________________________________________________
 955 Bool_t AliDielectronMC::CheckParticleSource(Int_t label, AliDielectronSignalMC::ESource source) const {
 956   //
 957   //  Check the source for the particle
 958   //  NOTE: for tracks the absolute label should be passed
 959   //
 960
 961   switch (source) {
 962     case AliDielectronSignalMC::kDontCare :
 963       return kTRUE;
 964     break;
 965     case AliDielectronSignalMC::kPrimary :
 966       // true if label is in the list of particles from physics generator
 967       // NOTE: This includes all physics event history (initial state particles,
 968       //       exchange bosons, quarks, di-quarks, strings, un-stable particles, final state particles)
 969       //       Only the final state particles make it to the detector!!
 970       return (label>=0 && label<=GetNPrimary());
 971     break;
 972     case AliDielectronSignalMC::kFinalState :
 973       // primary particles created in the collision which reach the detectors
 974       // These would be:
 975       // 1.) particles produced in the collision
 976       // 2.) stable particles with respect to strong and electromagnetic interactions
 977       // 3.) excludes initial state particles
 978       // 4.) includes products of directly produced Sigma0 hyperon decay
 979       // 5.) includes products of directly produced pi0 decays
 980       // 6.) includes products of directly produced beauty hadron decays
 981       return IsPhysicalPrimary(label);
 982     break;
 983     case AliDielectronSignalMC::kDirect :
 984       // Primary particles which do not have any mother
 985       // This is the case for:
 986       // 1.) Initial state particles (the 2 protons in Pythia pp collisions)
 987       // 2.) In some codes, with sudden freeze-out, all particles generated from the fireball are direct.
 988       //     There is no history for these particles.
 989       // 3.) Certain particles added via MC generator cocktails (e.g. J/psi added to pythia MB events)
 990       return (label>=0 && GetMothersLabel(label)<0);
 991       break;
 992     case AliDielectronSignalMC::kNoCocktail :
 993       // Particles from the HIJING event and NOT from the AliGenCocktail
 994       return (label>=0 && GetMothersLabel(label)>=0);
 995       break;
 996     case AliDielectronSignalMC::kSecondary :
 997       // particles which are created by the interaction of final state primaries with the detector
 998       // or particles from strange weakly decaying particles (e.g. lambda, kaons, etc.)
 999       return (label>=GetNPrimary() && !IsPhysicalPrimary(label));
1000     break;
1001     default :
1002       return kFALSE;
1003   }
1004   return kFALSE;
1005 }
1006
1007 //________________________________________________________________________________
1008 Bool_t AliDielectronMC::CheckIsRadiative(Int_t label) const
1009 {
1010   //
1011   // Check if the particle has a three body decay, one being a photon
1012   //
1013   if(label<0) return kFALSE;
1014
1015
1016   if(fAnaType==kAOD) {
1017     if(!fMcArray) return kFALSE;
1018     AliAODMCParticle *mother=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label));
1019     if (!mother) return kFALSE;
1020     const Int_t nd=mother->GetNDaughters();
1021     if (nd==2) return kFALSE;
1022     for (Int_t i=2; i<nd; ++i)
1023       if (GetMCTrackFromMCEvent(mother->GetDaughter(0)+i)->PdgCode()!=22) return kFALSE; //last daughter is photon
1024   } else if(fAnaType==kESD) {
1025     if (!fMCEvent) return kFALSE;
1026     AliMCParticle *mother=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(label));
1027     if (!mother) return kFALSE;
1028     const Int_t nd=(mother->GetLastDaughter()-mother->GetFirstDaughter()+1);
1029     if (nd==2) return kFALSE;
1030     for (Int_t i=2; i<nd; ++i)
1031       if (GetMCTrackFromMCEvent(mother->GetFirstDaughter()+i)->PdgCode()!=22) return kFALSE; //last daughters are photons
1032   }
1033   return kTRUE;
1034 }
1035
1036 //________________________________________________________________________________
1037 Bool_t AliDielectronMC::CheckRadiativeDecision(Int_t mLabel, const AliDielectronSignalMC * const signalMC) const
1038 {
1039   //
1040   // Check for the decision of the radiative type request
1041   //
1042
1043   if (!signalMC) return kFALSE;
1044
1045   if (signalMC->GetJpsiRadiative()==AliDielectronSignalMC::kAll) return kTRUE;
1046
1047   Bool_t isRadiative=CheckIsRadiative(mLabel);
1048   if ((signalMC->GetJpsiRadiative()==AliDielectronSignalMC::kIsRadiative) && !isRadiative) return kFALSE;
1049   if ((signalMC->GetJpsiRadiative()==AliDielectronSignalMC::kIsNotRadiative) && isRadiative) return kFALSE;
1050
1051   return kTRUE;
1052 }
1053
1054 //________________________________________________________________________________
1055 Bool_t AliDielectronMC::IsMCTruth(Int_t label, AliDielectronSignalMC* signalMC, Int_t branch) const {
1056   //
1057   // Check if the particle corresponds to the MC truth in signalMC in the branch specified
1058   //
1059
1060   // NOTE:  Some particles have the sign of the label flipped. It is related to the quality of matching
1061   //        between the ESD and the MC track. The negative labels indicate a poor matching quality
1062   //if(label<0) return kFALSE;
1063   if(label<0) label *= -1;
1064
1065   AliVParticle* part = GetMCTrackFromMCEvent(label);
1066   if (!part) {
1067     AliError(Form("Could not find MC particle with label %d",label));
1068     return kFALSE;
1069   }
1070
1071   // check the leg
1072   if(!ComparePDG(part->PdgCode(),signalMC->GetLegPDG(branch),signalMC->GetLegPDGexclude(branch),signalMC->GetCheckBothChargesLegs(branch))) return kFALSE;
1073   if(!CheckParticleSource(label, signalMC->GetLegSource(branch))) return kFALSE;
1074
1075   // check the mother
1076   AliVParticle* mcMother=0x0;
1077   Int_t mLabel = -1;
1078   if(signalMC->GetMotherPDG(branch)!=0 || signalMC->GetMotherSource(branch)!=AliDielectronSignalMC::kDontCare) {
1079     if(part) {
1080       mLabel = GetMothersLabel(label);
1081       mcMother = GetMCTrackFromMCEvent(mLabel);
1082     }
1083     if(!mcMother && !signalMC->GetMotherPDGexclude(branch)) return kFALSE;
1084
1085     if(!ComparePDG((mcMother ? mcMother->PdgCode() : 0),signalMC->GetMotherPDG(branch),signalMC->GetMotherPDGexclude(branch),signalMC->GetCheckBothChargesMothers(branch))) return kFALSE;
1086     if(!CheckParticleSource(mLabel, signalMC->GetMotherSource(branch))) return kFALSE;
1087
1088     //check for radiative deday
1089     if (!CheckRadiativeDecision(mLabel, signalMC)) return kFALSE;
1090   }
1091
1092   // check the grandmother
1093   AliVParticle* mcGrandMother=0x0;
1094   if(signalMC->GetGrandMotherPDG(branch)!=0 || signalMC->GetGrandMotherSource(branch)!=AliDielectronSignalMC::kDontCare) {
1095     Int_t gmLabel = -1;
1096     if(mcMother) {
1097       gmLabel = GetMothersLabel(mLabel);
1098       mcGrandMother = static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(gmLabel));
1099     }
1100     if(!mcGrandMother && !signalMC->GetGrandMotherPDGexclude(branch)) return kFALSE;
1101
1102     if(!ComparePDG((mcGrandMother ? mcGrandMother->PdgCode() : 0),signalMC->GetGrandMotherPDG(branch),signalMC->GetGrandMotherPDGexclude(branch),signalMC->GetCheckBothChargesGrandMothers(branch))) return kFALSE;
1103     if(!CheckParticleSource(gmLabel, signalMC->GetGrandMotherSource(branch))) return kFALSE;
1104   }
1105
1106   return kTRUE;
1107 }
1108
1109
1110
1111 //________________________________________________________________________________
1112 Bool_t AliDielectronMC::IsMCTruth(const AliDielectronPair* pair, const AliDielectronSignalMC* signalMC) const {
1113   //
1114   // Check if the pair corresponds to the MC truth in signalMC
1115   //
1116
1117   // legs (daughters)
1118   const AliVParticle * mcD1 = pair->GetFirstDaughterP();
1119   const AliVParticle * mcD2 = pair->GetSecondDaughterP();
1120   Int_t labelD1 = (mcD1 ? TMath::Abs(mcD1->GetLabel()) : -1);
1121   Int_t labelD2 = (mcD2 ? TMath::Abs(mcD2->GetLabel()) : -1);
1122   Int_t d1Pdg = GetPdgFromLabel(labelD1);
1123   Int_t d2Pdg = GetPdgFromLabel(labelD2);
1124
1125   // mothers
1126   AliVParticle* mcM1=0x0;
1127   AliVParticle* mcM2=0x0;
1128
1129   // grand-mothers
1130   AliVParticle* mcG1 = 0x0;
1131   AliVParticle* mcG2 = 0x0;
1132
1133   // make direct(1-1 and 2-2) and cross(1-2 and 2-1) comparisons for the whole branch
1134   Bool_t directTerm = kTRUE;
1135   // daughters
1136   directTerm = directTerm && mcD1 && ComparePDG(d1Pdg,signalMC->GetLegPDG(1),signalMC->GetLegPDGexclude(1),signalMC->GetCheckBothChargesLegs(1))
1137                && CheckParticleSource(labelD1, signalMC->GetLegSource(1));
1138
1139   directTerm = directTerm && mcD2 && ComparePDG(d2Pdg,signalMC->GetLegPDG(2),signalMC->GetLegPDGexclude(2),signalMC->GetCheckBothChargesLegs(2))
1140                && CheckParticleSource(labelD2, signalMC->GetLegSource(2));
1141
1142   // mothers
1143   Int_t labelM1 = -1;
1144   if(signalMC->GetMotherPDG(1)!=0 || signalMC->GetMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
1145     labelM1 = GetMothersLabel(labelD1);
1146     if(labelD1>-1 && labelM1>-1) mcM1 = GetMCTrackFromMCEvent(labelM1);
1147     directTerm = directTerm && (mcM1 || signalMC->GetMotherPDGexclude(1))
1148                  && ComparePDG((mcM1 ? mcM1->PdgCode() : 0),signalMC->GetMotherPDG(1),signalMC->GetMotherPDGexclude(1),signalMC->GetCheckBothChargesMothers(1))
1149                  && CheckParticleSource(labelM1, signalMC->GetMotherSource(1))
1150                  && CheckRadiativeDecision(labelM1,signalMC);
1151   }
1152
1153   Int_t labelM2 = -1;
1154   if(signalMC->GetMotherPDG(2)!=0 || signalMC->GetMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
1155     labelM2 = GetMothersLabel(labelD2);
1156     if(labelD2>-1 && labelM2>-1) mcM2 = GetMCTrackFromMCEvent(labelM2);
1157     directTerm = directTerm && (mcM2 || signalMC->GetMotherPDGexclude(2))
1158                  && ComparePDG((mcM2 ? mcM2->PdgCode() : 0),signalMC->GetMotherPDG(2),signalMC->GetMotherPDGexclude(2),signalMC->GetCheckBothChargesMothers(2))
1159                  && CheckParticleSource(labelM2, signalMC->GetMotherSource(2))
1160                  && CheckRadiativeDecision(labelM2,signalMC);
1161   }
1162
1163   // grand-mothers
1164   Int_t labelG1 = -1;
1165   if(signalMC->GetGrandMotherPDG(1)!=0 || signalMC->GetGrandMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
1166     labelG1 = GetMothersLabel(labelM1);
1167     if(mcM1 && labelG1>-1) mcG1 = GetMCTrackFromMCEvent(labelG1);
1168     directTerm = directTerm && (mcG1 || signalMC->GetGrandMotherPDGexclude(1))
1169                  && ComparePDG((mcG1 ? mcG1->PdgCode() : 0),signalMC->GetGrandMotherPDG(1),signalMC->GetGrandMotherPDGexclude(1),signalMC->GetCheckBothChargesGrandMothers(1))
1170                  && CheckParticleSource(labelG1, signalMC->GetGrandMotherSource(1));
1171   }
1172
1173   Int_t labelG2 = -1;
1174   if(signalMC->GetGrandMotherPDG(2)!=0 || signalMC->GetGrandMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
1175     labelG2 = GetMothersLabel(labelM2);
1176     if(mcM2 && labelG2>-1) mcG2 = GetMCTrackFromMCEvent(labelG2);
1177     directTerm = directTerm && (mcG2 || signalMC->GetGrandMotherPDGexclude(2))
1178                  && ComparePDG((mcG2 ? mcG2->PdgCode() : 0),signalMC->GetGrandMotherPDG(2),signalMC->GetGrandMotherPDGexclude(2),signalMC->GetCheckBothChargesGrandMothers(2))
1179                  && CheckParticleSource(labelG2, signalMC->GetGrandMotherSource(2));
1180   }
1181
1182   // Cross term
1183   Bool_t crossTerm = kTRUE;
1184   // daughters
1185   crossTerm = crossTerm && mcD2
1186               && ComparePDG(d2Pdg,signalMC->GetLegPDG(1),signalMC->GetLegPDGexclude(1),signalMC->GetCheckBothChargesLegs(1))
1187               && CheckParticleSource(labelD2, signalMC->GetLegSource(1));
1188
1189   crossTerm = crossTerm && mcD1
1190               && ComparePDG(d1Pdg,signalMC->GetLegPDG(2),signalMC->GetLegPDGexclude(2),signalMC->GetCheckBothChargesLegs(2))
1191               && CheckParticleSource(labelD1, signalMC->GetLegSource(2));
1192
1193   // mothers
1194   if(signalMC->GetMotherPDG(1)!=0 || signalMC->GetMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
1195     if(!mcM2 && labelD2>-1) {
1196       labelM2 = GetMothersLabel(labelD2);
1197       if(labelM2>-1) mcM2 = GetMCTrackFromMCEvent(labelM2);
1198     }
1199     crossTerm = crossTerm && (mcM2 || signalMC->GetMotherPDGexclude(1))
1200                 && ComparePDG((mcM2 ? mcM2->PdgCode() : 0),signalMC->GetMotherPDG(1),signalMC->GetMotherPDGexclude(1),signalMC->GetCheckBothChargesMothers(1))
1201                 && CheckParticleSource(labelM2, signalMC->GetMotherSource(1))
1202                 && CheckRadiativeDecision(labelM2,signalMC);
1203   }
1204
1205   if(signalMC->GetMotherPDG(2)!=0 || signalMC->GetMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
1206     if(!mcM1 && labelD1>-1) {
1207       labelM1 = GetMothersLabel(labelD1);
1208       if(labelM1>-1) mcM1 = GetMCTrackFromMCEvent(labelM1);
1209     }
1210     crossTerm = crossTerm && (mcM1 || signalMC->GetMotherPDGexclude(2))
1211                 && ComparePDG((mcM1 ? mcM1->PdgCode() : 0),signalMC->GetMotherPDG(2),signalMC->GetMotherPDGexclude(2),signalMC->GetCheckBothChargesMothers(2))
1212                 && CheckParticleSource(labelM1, signalMC->GetMotherSource(2))
1213                 && CheckRadiativeDecision(labelM1,signalMC);
1214   }
1215
1216   // grand-mothers
1217   if(signalMC->GetGrandMotherPDG(1)!=0 || signalMC->GetGrandMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
1218     if(!mcG2 && mcM2) {
1219       labelG2 = GetMothersLabel(labelM2);
1220       if(labelG2>-1) mcG2 = GetMCTrackFromMCEvent(labelG2);
1221     }
1222     crossTerm = crossTerm && (mcG2 || signalMC->GetGrandMotherPDGexclude(1))
1223                 && ComparePDG((mcG2 ? mcG2->PdgCode() : 0),signalMC->GetGrandMotherPDG(1),signalMC->GetGrandMotherPDGexclude(1),signalMC->GetCheckBothChargesGrandMothers(1))
1224                 && CheckParticleSource(labelG2, signalMC->GetGrandMotherSource(1));
1225   }
1226
1227   if(signalMC->GetGrandMotherPDG(2)!=0 || signalMC->GetGrandMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
1228     if(!mcG1 && mcM1) {
1229       labelG1 = GetMothersLabel(labelM1);
1230       if(labelG1>-1) mcG1 = GetMCTrackFromMCEvent(labelG1);
1231     }
1232     crossTerm = crossTerm && (mcG1 || signalMC->GetGrandMotherPDGexclude(2))
1233                 && ComparePDG((mcG1 ? mcG1->PdgCode() : 0),signalMC->GetGrandMotherPDG(2),signalMC->GetGrandMotherPDGexclude(2),signalMC->GetCheckBothChargesGrandMothers(2))
1234                 && CheckParticleSource(labelG1, signalMC->GetGrandMotherSource(2));
1235   }
1236
1237   Bool_t motherRelation = kTRUE;
1238   if(signalMC->GetMothersRelation()==AliDielectronSignalMC::kSame) {
1239     motherRelation = motherRelation && HaveSameMother(pair);
1240   }
1241   if(signalMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent) {
1242     motherRelation = motherRelation && !HaveSameMother(pair);
1243   }
1244
1245   return ((directTerm || crossTerm) && motherRelation);
1246 }
1247
1248
1249
1250 //____________________________________________________________
1251 Bool_t AliDielectronMC::HaveSameMother(const AliDielectronPair * pair) const
1252 {
1253   //
1254   // Check whether two particles have the same mother
1255   //
1256
1257   const AliVParticle * daughter1 = pair->GetFirstDaughterP();
1258   const AliVParticle * daughter2 = pair->GetSecondDaughterP();
1259   if (!daughter1 || !daughter2) return 0;
1260
1261   AliVParticle *mcDaughter1=GetMCTrackFromMCEvent(daughter1->GetLabel());
1262   AliVParticle *mcDaughter2=GetMCTrackFromMCEvent(daughter2->GetLabel());
1263   if (!mcDaughter1 || !mcDaughter2) return 0;
1264
1265   Int_t labelMother1=-1;
1266   Int_t labelMother2=-1;
1267
1268   if (mcDaughter1->IsA()==AliMCParticle::Class()){
1269     labelMother1=(static_cast<AliMCParticle*>(mcDaughter1))->GetMother();
1270     labelMother2=(static_cast<AliMCParticle*>(mcDaughter2))->GetMother();
1271   } else if (mcDaughter1->IsA()==AliAODMCParticle::Class()) {
1272     labelMother1=(static_cast<AliAODMCParticle*>(mcDaughter1))->GetMother();
1273     labelMother2=(static_cast<AliAODMCParticle*>(mcDaughter2))->GetMother();
1274   }
1275
1276   Bool_t sameMother=(labelMother1>-1)&&(labelMother2>-1)&&(labelMother1==labelMother2);
1277
1278   return sameMother;
1279 }
1280
1281 //________________________________________________________________
1282 Int_t AliDielectronMC::IsJpsiPrimary(const AliDielectronPair * pair)
1283 {
1284  // return: "0" for primary jpsi
1285  //         "1" for secondary jpsi (from beauty)
1286  //         "2" for background
1287  if(!HaveSameMother(pair)) return 2;
1288  AliVParticle *mcDaughter1=GetMCTrackFromMCEvent((pair->GetFirstDaughterP())->GetLabel());
1289  Int_t labelMother=-1;
1290
1291   if (mcDaughter1->IsA()==AliMCParticle::Class()){
1292      labelMother=(static_cast<AliMCParticle*>(mcDaughter1))->GetMother();
1293      } else if (mcDaughter1->IsA()==AliAODMCParticle::Class()) {
1294      labelMother=(static_cast<AliAODMCParticle*>(mcDaughter1))->GetMother();
1295      }
1296
1297  AliVParticle* mcMother=GetMCTrackFromMCEvent(labelMother);
1298  if(!IsMCMotherToEE(mcMother,443)) return 2;
1299  return IsJpsiPrimary(mcMother);
1300 }
1301
1302 //______________________________________________________________
1303 Int_t AliDielectronMC::IsJpsiPrimary(const AliVParticle * particle)
1304 {
1305   // return: "0" for primary jpsi
1306   //         "1" for secondary jpsi (come from B decay)
1307  Int_t labelMoth=-1;
1308  Int_t pdgCode;
1309
1310  if (particle->IsA()==AliMCParticle::Class()){
1311      labelMoth = (static_cast<const AliMCParticle*>(particle))->GetMother();
1312      while(labelMoth>0){
1313        particle = GetMCTrackFromMCEvent(labelMoth);
1314        pdgCode = TMath::Abs((static_cast<const AliMCParticle*>(particle))->PdgCode());
1315        if((pdgCode>500 && pdgCode<600) || (pdgCode>5000 && pdgCode<6000)) return 1;
1316        labelMoth = (static_cast<const AliMCParticle*>(particle))->GetMother();
1317        }
1318     }
1319  else if (particle->IsA()==AliAODMCParticle::Class()){
1320      labelMoth = (static_cast<const AliAODMCParticle*>(particle))->GetMother();
1321      while(labelMoth>0){
1322      particle = GetMCTrackFromMCEvent(labelMoth);
1323      pdgCode = TMath::Abs((static_cast<const AliAODMCParticle*>(particle))->PdgCode());
1324      if((pdgCode>500 && pdgCode<600) || (pdgCode>5000 && pdgCode<6000)) return 1;
1325      labelMoth = (static_cast<const AliAODMCParticle*>(particle))->GetMother();
1326      }
1327   }
1328   return 0;
1329 }
1330
1331
1332 Bool_t AliDielectronMC::GetPrimaryVertex(Double_t &primVtxX, Double_t &primVtxY, Double_t &primVtxZ){
1333
1334      if(fAnaType == kESD){
1335      const AliVVertex* mcVtx =  fMCEvent->GetPrimaryVertex();
1336      if(!mcVtx) return kFALSE;
1337      primVtxX = mcVtx->GetX();
1338      primVtxY = mcVtx->GetY();
1339      primVtxZ = mcVtx->GetZ();
1340      }else if(fAnaType == kAOD){
1341      AliAODEvent *aod=((AliAODInputHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()))->GetEvent();
1342      if(!aod) return kFALSE;
1343      AliAODMCHeader *mcHead = dynamic_cast<AliAODMCHeader*>(aod->FindListObject(AliAODMCHeader::StdBranchName()));
1344      if(!mcHead) return kFALSE;
1345      primVtxX = mcHead->GetVtxX();
1346      primVtxY = mcHead->GetVtxY();
1347      primVtxZ = mcHead->GetVtxZ();
1348      }
1349
1350 return kTRUE;
1351 }