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