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