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