PWGDQ/dielectron/AliDielectronMC.cxx

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