1 /*************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //#####################################################
18 //# Class AliDielectronMC #
19 //# Cut Class for Jpsi->e+e- analysis #
21 //# by WooJin J. Park, GSI / W.J.Park@gsi.de #
23 //#####################################################
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>
35 #include <AliESDEvent.h>
36 #include <AliAODEvent.h>
37 #include <AliESDtrack.h>
38 #include <AliAODTrack.h>
41 #include <TClonesArray.h>
42 #include <TParticle.h>
44 #include "AliDielectronSignalMC.h"
45 #include "AliDielectronMC.h"
48 ClassImp(AliDielectronMC)
50 AliDielectronMC* AliDielectronMC::fgInstance=0x0;
52 //____________________________________________________________
53 AliDielectronMC* AliDielectronMC::Instance()
56 // return pointer to singleton implementation
58 if (fgInstance) return fgInstance;
60 AnalysisType type=kUNSET;
62 if (AliAnalysisManager::GetAnalysisManager()){
63 if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliESDInputHandler::Class()) type=kESD;
64 else if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliAODInputHandler::Class()) type=kAOD;
66 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
67 if(type == kESD) hasMC=mcHandler!=0x0;
70 fgInstance=new AliDielectronMC(type);
72 fgInstance->SetHasMC(hasMC);
77 //____________________________________________________________
78 AliDielectronMC::AliDielectronMC(AnalysisType type):
86 // default constructor
91 //____________________________________________________________
92 AliDielectronMC::~AliDielectronMC()
100 //____________________________________________________________
101 void AliDielectronMC::Initialize()
104 // initialize MC class
106 if (!ConnectMCEvent()) AliError("Initialization of MC object failed!");
109 //____________________________________________________________
110 Int_t AliDielectronMC::GetNMCTracks()
113 // return the number of generated tracks from MC event
115 if(fAnaType == kESD){
116 if (!fMCEvent){ AliError("No fMCEvent"); return 0; }
117 return fMCEvent->GetNumberOfTracks();}
118 else if(fAnaType == kAOD){
119 if(!fMcArray) { AliError("No fMcArray"); return 0; }
120 return fMcArray->GetEntriesFast();
125 //____________________________________________________________
126 Int_t AliDielectronMC::GetNMCTracksFromStack()
129 // return the number of generated tracks from stack
131 if (!fStack){ AliError("No fStack"); return -999; }
132 return fStack->GetNtrack();
135 //____________________________________________________________
136 Int_t AliDielectronMC::GetNPrimary() const
139 // return the number of primary track from MC event
141 if (!fMCEvent){ AliError("No fMCEvent"); return 0; }
142 return fMCEvent->GetNumberOfPrimaries();
145 //____________________________________________________________
146 Int_t AliDielectronMC::GetNPrimaryFromStack()
149 // return the number of primary track from stack
151 if (!fStack){ AliError("No fStack"); return -999; }
152 return fStack->GetNprimary();
155 //____________________________________________________________
156 AliVParticle* AliDielectronMC::GetMCTrackFromMCEvent(Int_t label) const
159 // return MC track directly from MC event
160 // used not only for tracks but for mothers as well, therefore do not use abs(label)
162 if (label<0) return NULL;
163 AliVParticle * track=0x0;
164 if(fAnaType == kESD){
165 if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
166 track = fMCEvent->GetTrack(label); // tracks from MC event (ESD)
167 } else if(fAnaType == kAOD) {
168 if (!fMcArray){ AliError("No fMcArray"); return NULL;}
169 if (label>fMcArray->GetEntriesFast()) { AliDebug(10,Form("track %d out of array size %d",label,fMcArray->GetEntriesFast())); return NULL;}
170 track = (AliVParticle*)fMcArray->At(label); // tracks from MC event (AOD)
175 //____________________________________________________________
176 Bool_t AliDielectronMC::ConnectMCEvent()
179 // connect stack object from the mc handler
185 if(fAnaType == kESD){
186 AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
187 if (!mcHandler){ /*AliError("Could not retrive MC event handler!");*/ return kFALSE; }
188 if (!mcHandler->InitOk() ) return kFALSE;
189 if (!mcHandler->TreeK() ) return kFALSE;
190 if (!mcHandler->TreeTR() ) return kFALSE;
192 AliMCEvent* mcEvent = mcHandler->MCEvent();
193 if (!mcEvent){ /*AliError("Could not retrieve MC event!");*/ return kFALSE; }
196 if (!UpdateStack()) return kFALSE;
198 else if(fAnaType == kAOD)
200 AliAODInputHandler* aodHandler=(AliAODInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
201 if (!aodHandler) return kFALSE;
202 AliAODEvent *aod=aodHandler->GetEvent();
203 if (!aod) return kFALSE;
205 fMcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
206 if (!fMcArray){ /*AliError("Could not retrieve MC array!");*/ return kFALSE; }
212 //____________________________________________________________
213 Bool_t AliDielectronMC::UpdateStack()
216 // update stack with new event
218 if (!fMCEvent){ AliError("No fMCEvent"); return kFALSE;}
219 AliStack* stack = fMCEvent->Stack();
220 if (!stack){ AliError("Could not retrive stack!"); return kFALSE; }
225 //____________________________________________________________
226 AliMCParticle* AliDielectronMC::GetMCTrack( const AliESDtrack* _track)
231 if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
232 Int_t nStack = fMCEvent->GetNumberOfTracks();
233 Int_t label = TMath::Abs(_track->GetLabel()); // negative label indicate poor matching quality
234 if(label>nStack)return NULL;
235 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(label));
240 //____________________________________________________________
241 AliAODMCParticle* AliDielectronMC::GetMCTrack( const AliAODTrack* _track)
246 if(!fMcArray) { AliError("No fMcArray"); return NULL;}
247 Int_t nStack = fMcArray->GetEntriesFast();
248 Int_t label = TMath::Abs(_track->GetLabel()); // negative label indicate poor matching quality
249 if(label > nStack) return NULL;
250 AliAODMCParticle *mctrack = (AliAODMCParticle*)fMcArray->At(label);
254 //____________________________________________________________
255 TParticle* AliDielectronMC::GetMCTrackFromStack(const AliESDtrack* _track)
258 // return MC track from stack
260 Int_t label = TMath::Abs(_track->GetLabel());
261 if (!fStack) AliWarning("fStack is not available. Update stack first.");
262 TParticle* mcpart = fStack->Particle(label);
263 if (!mcpart) return NULL;
267 //____________________________________________________________
268 AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliESDtrack* _track)
271 // return MC track mother
273 AliMCParticle* mcpart = GetMCTrack(_track);
274 if (!mcpart) return NULL;
275 if(mcpart->GetMother()<0) return NULL;
276 AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(mcpart->GetMother()));
277 if (!mcmother) return NULL;
281 //______________________________________________________________
282 AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODTrack* _track)
285 // return MC track mother
287 AliAODMCParticle* mcpart = GetMCTrack(_track);
288 if (!mcpart) return NULL;
289 if(mcpart->GetMother() < 0) return NULL;
290 AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMcArray->At(mcpart->GetMother()));
291 if (!mcmother) return NULL;
294 //____________________________________________________________
295 AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliMCParticle* _particle){
297 // return MC track mother
299 if(_particle->GetMother() < 0) return NULL;
300 AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(_particle->GetMother()));
304 //____________________________________________________________
305 AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODMCParticle* _particle){
307 // return MC track mother
309 if( _particle->GetMother() < 0) return NULL;
310 AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMcArray->At(_particle->GetMother()));
314 //____________________________________________________________
315 TParticle* AliDielectronMC::GetMCTrackMotherFromStack(const AliESDtrack* _track)
318 // return MC track mother from stack
320 TParticle* mcpart = GetMCTrackFromStack(_track);
321 if ( !mcpart || mcpart->GetFirstMother()<=0 ) return NULL;
322 TParticle* mcmother = fStack->Particle(mcpart->GetFirstMother());
323 if (!mcmother) return NULL;
327 //____________________________________________________________
328 Int_t AliDielectronMC::GetMCPID(const AliESDtrack* _track)
331 // return PDG code of the track from the MC truth info
333 AliMCParticle* mcpart = GetMCTrack(_track);
334 if (!mcpart) return -999;
335 return mcpart->PdgCode();
338 //__________________________________________________________
339 Int_t AliDielectronMC::GetMCPID(const AliAODTrack* _track)
342 // return PDG code of the track from the MC truth info
344 AliAODMCParticle* mcpart = GetMCTrack(_track);
345 if (!mcpart) return -999;
346 return mcpart->PdgCode();
349 //____________________________________________________________
350 Int_t AliDielectronMC::GetMCPIDFromStack(const AliESDtrack* _track)
353 // return MC PDG code from stack
355 TParticle* mcpart = GetMCTrackFromStack(_track);
356 if (!mcpart) return -999;
357 return mcpart->GetPdgCode();
360 //____________________________________________________________
361 Int_t AliDielectronMC::GetMotherPDG( const AliESDtrack* _track)
364 // return PDG code of the mother track from the MC truth info
366 AliMCParticle* mcmother = GetMCTrackMother(_track);
367 if (!mcmother) return -999;
368 return mcmother->PdgCode();
371 //________________________________________________________
372 Int_t AliDielectronMC::GetMotherPDG( const AliAODTrack* _track)
375 // return PDG code of the mother track from the MC truth info
377 AliAODMCParticle* mcmother = GetMCTrackMother(_track);
378 if (!mcmother) return -999;
379 return mcmother->PdgCode();
382 //________________________________________________________
383 Int_t AliDielectronMC::GetMotherPDG( const AliMCParticle* _track)
386 // return PDG code of the mother track from the MC truth info
388 AliMCParticle* mcmother = GetMCTrackMother(_track);
389 if (!mcmother) return -999;
390 return mcmother->PdgCode();
393 //________________________________________________________
394 Int_t AliDielectronMC::GetMotherPDG( const AliAODMCParticle* _track)
397 // return PDG code of the mother track from the MC truth info
399 AliAODMCParticle* mcmother = GetMCTrackMother(_track);
400 if (!mcmother) return -999;
401 return mcmother->PdgCode();
404 //____________________________________________________________
405 Int_t AliDielectronMC::GetMotherPDGFromStack(const AliESDtrack* _track)
408 // return PDG code of the mother track from stack
410 TParticle* mcmother = GetMCTrackMotherFromStack(_track);
411 if (!mcmother) return -999;
412 return mcmother->GetPdgCode();
415 //____________________________________________________________
416 Int_t AliDielectronMC::GetMCProcess(const AliESDtrack* _track)
419 // return process number of the track
421 AliMCParticle* mcpart = GetMCTrack(_track);
422 if (!mcpart) return -999;
426 //____________________________________________________________
427 Int_t AliDielectronMC::GetMCProcessFromStack(const AliESDtrack* _track)
430 // return process number of the track
432 TParticle* mcpart = GetMCTrackFromStack(_track);
433 if (!mcpart) return -999;
434 return mcpart->GetUniqueID();
437 //____________________________________________________________
438 Int_t AliDielectronMC::NumberOfDaughters(const AliESDtrack* track)
441 // returns the number of daughters
443 AliMCParticle *mcmother=GetMCTrackMother(track);
444 if(!mcmother||!mcmother->Particle()) return -999;
445 // return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
446 return mcmother->Particle()->GetNDaughters();
449 //_________________________________________________________
450 Int_t AliDielectronMC::NumberOfDaughters(const AliAODTrack* track)
453 // returns the number of daughters
455 AliAODMCParticle *mcmother=GetMCTrackMother(track);
456 if(!mcmother) return -999;
457 return NumberOfDaughters(mcmother);
461 //____________________________________________________________
462 Int_t AliDielectronMC::NumberOfDaughters(const AliMCParticle* particle)
465 // returns the number of daughters
467 AliMCParticle *mcmother=GetMCTrackMother(particle);
468 if(!mcmother||!mcmother->Particle()) return -999;
469 //return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
470 return mcmother->Particle()->GetNDaughters();
473 //____________________________________________________________
474 Int_t AliDielectronMC::NumberOfDaughters(const AliAODMCParticle* particle)
477 // returns the number of daughters
479 AliAODMCParticle *mcmother=GetMCTrackMother(particle);
480 if(!mcmother) return -999;
481 return mcmother->GetNDaughters();
484 //____________________________________________________________
485 Int_t AliDielectronMC::GetMCProcessMother(const AliESDtrack* _track)
488 // return process number of the mother of the track
490 AliMCParticle* mcmother = GetMCTrackMother(_track);
491 if (!mcmother) return -999;
495 //____________________________________________________________
496 Int_t AliDielectronMC::GetMCProcessMotherFromStack(const AliESDtrack* _track)
499 // return process number of the mother of the track
501 TParticle* mcmother = GetMCTrackMotherFromStack(_track);
502 if (!mcmother) return -999;
503 return mcmother->GetUniqueID();
506 //____________________________________________________________
507 Bool_t AliDielectronMC::IsMCMotherToEE(const AliVParticle *particle, Int_t pdgMother)
510 // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
512 if (fAnaType==kESD && !fMCEvent) return kFALSE;
513 if (fAnaType==kAOD && !fMcArray) return kFALSE;
514 if (!particle) return kFALSE;
516 if (particle->IsA()==AliMCParticle::Class()){
517 return IsMCMotherToEEesd(static_cast<const AliMCParticle*>(particle),pdgMother);
518 } else if (particle->IsA()==AliAODMCParticle::Class()){
519 return IsMCMotherToEEaod(static_cast<const AliAODMCParticle*>(particle),pdgMother);
521 AliError("Unknown particle type");
526 //____________________________________________________________
527 Bool_t AliDielectronMC::IsMCMotherToEEesd(const AliMCParticle *particle, Int_t pdgMother)
530 // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
535 if (particle->PdgCode()!=pdgMother) return kFALSE;
536 Int_t ifirst = particle->GetFirstDaughter();
537 Int_t ilast = particle->GetLastDaughter();
539 //check number of daughters
540 if ((ilast-ifirst)!=1) return kFALSE;
541 AliMCParticle *firstD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ifirst));
542 AliMCParticle *secondD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ilast));
544 //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
545 if (firstD->Charge()>0){
546 if (firstD->PdgCode()!=-11) return kFALSE;
547 if (secondD->PdgCode()!=11) return kFALSE;
549 if (firstD->PdgCode()!=11) return kFALSE;
550 if (secondD->PdgCode()!=-11) return kFALSE;
556 //____________________________________________________________
557 Bool_t AliDielectronMC::IsMCMotherToEEaod(const AliAODMCParticle *particle, Int_t pdgMother)
560 // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
564 if (particle->GetPdgCode()!=pdgMother) return kFALSE;
565 if (particle->GetNDaughters()!=2) return kFALSE;
567 Int_t ifirst = particle->GetDaughter(0);
568 Int_t ilast = particle->GetDaughter(1);
570 //check number of daughters
571 if ((ilast-ifirst)!=1) return kFALSE;
573 AliAODMCParticle *firstD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ifirst));
574 AliAODMCParticle *secondD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ilast));
576 //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
578 if (firstD->Charge()>0){
579 if (firstD->GetPdgCode()!=-11) return kFALSE;
580 if (secondD->GetPdgCode()!=11) return kFALSE;
582 if (firstD->GetPdgCode()!=11) return kFALSE;
583 if (secondD->GetPdgCode()!=-11) return kFALSE;
588 //____________________________________________________________
589 Int_t AliDielectronMC::GetLabelMotherWithPdg(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
592 // test if mother of particle 1 and 2 has pdgCode pdgMother and is the same;
595 if (!fMCEvent) return -1;
596 return GetLabelMotherWithPdgESD(particle1, particle2, pdgMother);
598 else if (fAnaType==kAOD)
600 if (!fMcArray) return -1;
601 return GetLabelMotherWithPdgAOD(particle1, particle2, pdgMother);
607 //____________________________________________________________
608 Int_t AliDielectronMC::GetLabelMotherWithPdgESD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
611 // test if mother of particle 1 and 2 has pdgCode +-11 (electron),
612 // have the same mother and the mother had pdg code pdgMother
614 //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
616 // negative label indicate poor matching quality
617 Int_t lblPart1 = TMath::Abs(particle1->GetLabel());
618 Int_t lblPart2 = TMath::Abs(particle2->GetLabel());
619 AliMCParticle *mcPart1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(lblPart1));
620 AliMCParticle *mcPart2=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(lblPart2));
622 if (!mcPart1||!mcPart2) return -1;
624 Int_t lblMother1=mcPart1->GetMother();
625 Int_t lblMother2=mcPart2->GetMother();
626 AliMCParticle *mcMother1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(lblMother1));
628 if (!mcMother1) return -1;
629 if (lblMother1!=lblMother2) return -1;
630 if (TMath::Abs(mcPart1->PdgCode())!=11) return -1;
631 if (mcPart1->PdgCode()!=-mcPart2->PdgCode()) return -1;
632 if (mcMother1->PdgCode()!=pdgMother) return -1;
637 //____________________________________________________________
638 Int_t AliDielectronMC::GetLabelMotherWithPdgAOD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
641 // test if mother of particle 1 and 2 has pdgCode +-11 (electron),
642 // have the same mother and the mother had pdg code pdgMother
644 //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
646 // negative label indicate poor matching quality
647 Int_t lblPart1 = TMath::Abs(particle1->GetLabel());
648 Int_t lblPart2 = TMath::Abs(particle2->GetLabel());
649 AliAODMCParticle *mcPart1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(lblPart1));
650 AliAODMCParticle *mcPart2=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(lblPart2));
652 if (!mcPart1||!mcPart2) return -1;
654 Int_t lblMother1=mcPart1->GetMother();
655 Int_t lblMother2=mcPart2->GetMother();
656 AliAODMCParticle *mcMother1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(lblMother1));
658 if (!mcMother1) return -1;
659 if (lblMother1!=lblMother2) return -1;
660 if (TMath::Abs(mcPart1->GetPdgCode())!=11) return -1;
661 if (mcPart1->GetPdgCode()!=-mcPart2->GetPdgCode()) return -1;
662 if (mcMother1->GetPdgCode()!=pdgMother) return -1;
667 //____________________________________________________________
668 void AliDielectronMC::GetDaughters(const TObject *mother, AliVParticle* &d1, AliVParticle* &d2)
671 // Get First two daughters of the mother
678 if(!fMcArray) return;
679 const AliAODMCParticle *aodMother=static_cast<const AliAODMCParticle*>(mother);
680 lblD1=aodMother->GetDaughter(0);
681 lblD2=aodMother->GetDaughter(1);
682 d1 = (AliVParticle*)fMcArray->At(lblD1);
683 d2 = (AliVParticle*)fMcArray->At(lblD2);
684 } else if (fAnaType==kESD){
685 if (!fMCEvent) return;
686 const AliMCParticle *aodMother=static_cast<const AliMCParticle*>(mother);
687 lblD1=aodMother->GetFirstDaughter();
688 lblD2=aodMother->GetLastDaughter();
689 d1=fMCEvent->GetTrack(lblD1);
690 d2=fMCEvent->GetTrack(lblD2);
695 //________________________________________________________________________________
696 Int_t AliDielectronMC::GetMothersLabel(Int_t daughterLabel) const {
698 // Get the label of the mother for particle with label daughterLabel
699 // NOTE: for tracks, the absolute label should be passed
701 if(daughterLabel<0) return -1;
702 if (fAnaType==kAOD) {
703 if(!fMcArray) return -1;
704 if(GetMCTrackFromMCEvent(daughterLabel))
705 return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(daughterLabel)))->GetMother();
706 } else if(fAnaType==kESD) {
707 if (!fMCEvent) return -1;
708 if(GetMCTrackFromMCEvent(daughterLabel))
709 return (static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(daughterLabel)))->GetMother();
715 //________________________________________________________________________________
716 Int_t AliDielectronMC::GetPdgFromLabel(Int_t label) const {
718 // Get particle code using the label from stack
719 // NOTE: for tracks, the absolute label should be passed
721 if(label<0) return 0;
723 if(!fMcArray) return 0;
724 return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label)))->PdgCode();
725 } else if(fAnaType==kESD) {
726 if (!fMCEvent) return 0;
727 return (static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(label)))->PdgCode();
733 //________________________________________________________________________________
734 Bool_t AliDielectronMC::ComparePDG(Int_t particlePDG, Int_t requiredPDG, Bool_t pdgExclusion, Bool_t checkBothCharges) const {
736 // Test the PDG codes of particles with the required ones
738 Bool_t result = kTRUE;
739 Int_t absRequiredPDG = TMath::Abs(requiredPDG);
741 switch(absRequiredPDG) {
743 result = kTRUE; // PDG not required (any code will do fine)
745 case 100: // light flavoured mesons
747 result = TMath::Abs(particlePDG)>=100 && TMath::Abs(particlePDG)<=199;
749 if(requiredPDG>0) result = particlePDG>=100 && particlePDG<=199;
750 if(requiredPDG<0) result = particlePDG>=-199 && particlePDG<=-100;
753 case 1000: // light flavoured baryons
755 result = TMath::Abs(particlePDG)>=1000 && TMath::Abs(particlePDG)<=1999;
757 if(requiredPDG>0) result = particlePDG>=1000 && particlePDG<=1999;
758 if(requiredPDG<0) result = particlePDG>=-1999 && particlePDG<=-1000;
761 case 200: // light flavoured mesons
763 result = TMath::Abs(particlePDG)>=200 && TMath::Abs(particlePDG)<=299;
765 if(requiredPDG>0)result = particlePDG>=200 && particlePDG<=299;
766 if(requiredPDG<0)result = particlePDG>=-299 && particlePDG<=-200;
769 case 2000: // light flavoured baryons
771 result = TMath::Abs(particlePDG)>=2000 && TMath::Abs(particlePDG)<=2999;
773 if(requiredPDG>0) result = particlePDG>=2000 && particlePDG<=2999;
774 if(requiredPDG<0) result = particlePDG>=-2999 && particlePDG<=-2000;
777 case 300: // all strange mesons
779 result = TMath::Abs(particlePDG)>=300 && TMath::Abs(particlePDG)<=399;
781 if(requiredPDG>0) result = particlePDG>=300 && particlePDG<=399;
782 if(requiredPDG<0) result = particlePDG>=-399 && particlePDG<=-300;
785 case 3000: // all strange baryons
787 result = TMath::Abs(particlePDG)>=3000 && TMath::Abs(particlePDG)<=3999;
789 if(requiredPDG>0) result = particlePDG>=3000 && particlePDG<=3999;
790 if(requiredPDG<0) result = particlePDG>=-3999 && particlePDG<=-3000;
793 case 400: // all charmed mesons
795 result = TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=499;
797 if(requiredPDG>0) result = particlePDG>=400 && particlePDG<=499;
798 if(requiredPDG<0) result = particlePDG>=-499 && particlePDG<=-400;
801 case 401: // open charm mesons
803 result = TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=439;
805 if(requiredPDG>0) result = particlePDG>=400 && particlePDG<=439;
806 if(requiredPDG<0) result = particlePDG>=-439 && particlePDG<=-400;
809 case 402: // open charm mesons and baryons together
811 result = (TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=439) ||
812 (TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4399);
814 if(requiredPDG>0) result = (particlePDG>=400 && particlePDG<=439) ||
815 (particlePDG>=4000 && particlePDG<=4399);
816 if(requiredPDG<0) result = (particlePDG>=-439 && particlePDG<=-400) ||
817 (particlePDG>=-4399 && particlePDG<=-4000);
820 case 403: // all charm hadrons
822 result = (TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=499) ||
823 (TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4999);
825 if(requiredPDG>0) result = (particlePDG>=400 && particlePDG<=499) ||
826 (particlePDG>=4000 && particlePDG<=4999);
827 if(requiredPDG<0) result = (particlePDG>=-499 && particlePDG<=-400) ||
828 (particlePDG>=-4999 && particlePDG<=-4000);
831 case 4000: // all charmed baryons
833 result = TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4999;
835 if(requiredPDG>0) result = particlePDG>=4000 && particlePDG<=4999;
836 if(requiredPDG<0) result = particlePDG>=-4999 && particlePDG<=-4000;
839 case 4001: // open charm baryons
841 result = TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4399;
843 if(requiredPDG>0) result = particlePDG>=4000 && particlePDG<=4399;
844 if(requiredPDG<0) result = particlePDG>=-4399 && particlePDG<=-4000;
847 case 500: // all beauty mesons
849 result = TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=599;
851 if(requiredPDG>0) result = particlePDG>=500 && particlePDG<=599;
852 if(requiredPDG<0) result = particlePDG>=-599 && particlePDG<=-500;
855 case 501: // open beauty mesons
857 result = TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=549;
859 if(requiredPDG>0) result = particlePDG>=500 && particlePDG<=549;
860 if(requiredPDG<0) result = particlePDG>=-549 && particlePDG<=-500;
863 case 502: // open beauty mesons and baryons
865 result = (TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=549) ||
866 (TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5499);
868 if(requiredPDG>0) result = (particlePDG>=500 && particlePDG<=549) ||
869 (particlePDG>=5000 && particlePDG<=5499);
870 if(requiredPDG<0) result = (particlePDG>=-549 && particlePDG<=-500) ||
871 (particlePDG>=-5499 && particlePDG<=-5000);
874 case 503: // all beauty hadrons
876 result = (TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=599) ||
877 (TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5999);
879 if(requiredPDG>0) result = (particlePDG>=500 && particlePDG<=599) ||
880 (particlePDG>=5000 && particlePDG<=5999);
881 if(requiredPDG<0) result = (particlePDG>=-599 && particlePDG<=-500) ||
882 (particlePDG>=-5999 && particlePDG<=-5000);
885 case 5000: // all beauty baryons
887 result = TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5999;
889 if(requiredPDG>0) result = particlePDG>=5000 && particlePDG<=5999;
890 if(requiredPDG<0) result = particlePDG>=-5999 && particlePDG<=-5000;
893 case 5001: // open beauty baryons
895 result = TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5499;
897 if(requiredPDG>0) result = particlePDG>=5000 && particlePDG<=5499;
898 if(requiredPDG<0) result = particlePDG>=-5499 && particlePDG<=-5000;
901 case 902: // // open charm,beauty mesons and baryons together
903 result = (TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=439) ||
904 (TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4399) ||
905 (TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=549) ||
906 (TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5499);
908 if(requiredPDG>0) result = (particlePDG>=400 && particlePDG<=439) ||
909 (particlePDG>=4000 && particlePDG<=4399) ||
910 (particlePDG>=500 && particlePDG<=549) ||
911 (particlePDG>=5000 && particlePDG<=5499);
912 if(requiredPDG<0) result = (particlePDG>=-439 && particlePDG<=-400) ||
913 (particlePDG>=-4399 && particlePDG<=-4000) ||
914 (particlePDG>=-549 && particlePDG<=-500) ||
915 (particlePDG>=-5499 && particlePDG<=-5000);
918 default: // all specific cases
920 result = (absRequiredPDG==TMath::Abs(particlePDG));
922 result = (requiredPDG==particlePDG);
925 if(absRequiredPDG!=0 && pdgExclusion) result = !result;
929 //________________________________________________________________________________
930 Bool_t AliDielectronMC::IsPhysicalPrimary(Int_t label) const {
932 // Check if the particle with label "label" is a physical primary according to the
933 // definition in AliStack::IsPhysicalPrimary(Int_t label)
934 // Convention for being physical primary:
935 // 1.) particles produced in the collision
936 // 2.) stable particles with respect to strong and electromagnetic interactions
937 // 3.) excludes initial state particles
938 // 4.) includes products of directly produced Sigma0 hyperon decay
939 // 5.) includes products of directly produced pi0 decays
940 // 6.) includes products of directly produced beauty hadron decays
942 if(label<0) return kFALSE;
944 if(!fMcArray) return kFALSE;
945 return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label)))->IsPhysicalPrimary();
946 } else if(fAnaType==kESD) {
947 if (!fMCEvent) return kFALSE;
948 return fStack->IsPhysicalPrimary(label);
953 //________________________________________________________________________________
954 Bool_t AliDielectronMC::IsSecondaryFromWeakDecay(Int_t label) const {
956 // Check if the particle with label "label" is a physical secondary from weak decay according to the
957 // definition in AliStack::IsSecondaryFromWeakDecay(Int_t label)
959 if(label<0) return kFALSE;
961 if(!fMcArray) return kFALSE;
962 return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label)))->IsSecondaryFromWeakDecay();
963 } else if(fAnaType==kESD) {
964 if (!fMCEvent) return kFALSE;
965 return fStack->IsSecondaryFromWeakDecay(label);
970 //________________________________________________________________________________
971 Bool_t AliDielectronMC::IsSecondaryFromMaterial(Int_t label) const {
973 // Check if the particle with label "label" is a physical secondary from weak decay according to the
974 // definition in AliStack::IsSecondaryFromMaterial(Int_t label)
976 if(label<0) return kFALSE;
978 if(!fMcArray) return kFALSE;
979 return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label)))->IsSecondaryFromMaterial();
980 } else if(fAnaType==kESD) {
981 if (!fMCEvent) return kFALSE;
982 return fStack->IsSecondaryFromMaterial(label);
988 //________________________________________________________________________________
989 Bool_t AliDielectronMC::CheckParticleSource(Int_t label, AliDielectronSignalMC::ESource source) const {
991 // Check the source for the particle
992 // NOTE: for tracks the absolute label should be passed
996 case AliDielectronSignalMC::kDontCare :
999 case AliDielectronSignalMC::kPrimary :
1000 // true if label is in the list of particles from physics generator
1001 // NOTE: This includes all physics event history (initial state particles,
1002 // exchange bosons, quarks, di-quarks, strings, un-stable particles, final state particles)
1003 // Only the final state particles make it to the detector!!
1004 return (label>=0 && label<=GetNPrimary());
1006 case AliDielectronSignalMC::kFinalState :
1007 // primary particles created in the collision which reach the detectors
1009 // 1.) particles produced in the collision
1010 // 2.) stable particles with respect to strong and electromagnetic interactions
1011 // 3.) excludes initial state particles
1012 // 4.) includes products of directly produced Sigma0 hyperon decay
1013 // 5.) includes products of directly produced pi0 decays
1014 // 6.) includes products of directly produced beauty hadron decays
1015 return IsPhysicalPrimary(label);
1017 case AliDielectronSignalMC::kDirect :
1018 // Primary particles which do not have any mother
1019 // This is the case for:
1020 // 1.) Initial state particles (the 2 protons in Pythia pp collisions)
1021 // 2.) In some codes, with sudden freeze-out, all particles generated from the fireball are direct.
1022 // There is no history for these particles.
1023 // 3.) Certain particles added via MC generator cocktails (e.g. J/psi added to pythia MB events)
1024 return (label>=0 && GetMothersLabel(label)<0);
1026 case AliDielectronSignalMC::kNoCocktail :
1027 // Particles from the HIJING event and NOT from the AliGenCocktail
1028 return (label>=0 && GetMothersLabel(label)>=0);
1030 case AliDielectronSignalMC::kSecondary :
1031 // particles which are created by the interaction of final state primaries with the detector
1032 // or particles from strange weakly decaying particles (e.g. lambda, kaons, etc.)
1033 return (label>=GetNPrimary() && !IsPhysicalPrimary(label));
1035 case AliDielectronSignalMC::kSecondaryFromWeakDecay :
1036 // secondary particle from weak decay
1037 // or particles from strange weakly decaying particles (e.g. lambda, kaons, etc.)
1038 return (IsSecondaryFromWeakDecay(label));
1040 case AliDielectronSignalMC::kSecondaryFromMaterial :
1041 // secondary particle from material
1042 return (IsSecondaryFromMaterial(label));
1050 //________________________________________________________________________________
1051 Bool_t AliDielectronMC::CheckIsRadiative(Int_t label) const
1054 // Check if the particle has a three body decay, one being a photon
1056 if(label<0) return kFALSE;
1059 if(fAnaType==kAOD) {
1060 if(!fMcArray) return kFALSE;
1061 AliAODMCParticle *mother=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label));
1062 if (!mother) return kFALSE;
1063 const Int_t nd=mother->GetNDaughters();
1064 if (nd==2) return kFALSE;
1065 for (Int_t i=2; i<nd; ++i)
1066 if (GetMCTrackFromMCEvent(mother->GetDaughter(0)+i)->PdgCode()!=22) return kFALSE; //last daughter is photon
1067 } else if(fAnaType==kESD) {
1068 if (!fMCEvent) return kFALSE;
1069 AliMCParticle *mother=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(label));
1070 if (!mother) return kFALSE;
1071 const Int_t nd=(mother->GetLastDaughter()-mother->GetFirstDaughter()+1);
1072 if (nd==2) return kFALSE;
1073 for (Int_t i=2; i<nd; ++i)
1074 if (GetMCTrackFromMCEvent(mother->GetFirstDaughter()+i)->PdgCode()!=22) return kFALSE; //last daughters are photons
1079 //________________________________________________________________________________
1080 Bool_t AliDielectronMC::CheckRadiativeDecision(Int_t mLabel, const AliDielectronSignalMC * const signalMC) const
1083 // Check for the decision of the radiative type request
1086 if (!signalMC) return kFALSE;
1088 if (signalMC->GetJpsiRadiative()==AliDielectronSignalMC::kAll) return kTRUE;
1090 Bool_t isRadiative=CheckIsRadiative(mLabel);
1091 if ((signalMC->GetJpsiRadiative()==AliDielectronSignalMC::kIsRadiative) && !isRadiative) return kFALSE;
1092 if ((signalMC->GetJpsiRadiative()==AliDielectronSignalMC::kIsNotRadiative) && isRadiative) return kFALSE;
1097 //________________________________________________________________________________
1098 Bool_t AliDielectronMC::IsMCTruth(Int_t label, AliDielectronSignalMC* signalMC, Int_t branch) const {
1100 // Check if the particle corresponds to the MC truth in signalMC in the branch specified
1103 // NOTE: Some particles have the sign of the label flipped. It is related to the quality of matching
1104 // between the ESD and the MC track. The negative labels indicate a poor matching quality
1105 //if(label<0) return kFALSE;
1106 if(label<0) label *= -1;
1108 AliVParticle* part = GetMCTrackFromMCEvent(label);
1110 AliError(Form("Could not find MC particle with label %d",label));
1115 if(!ComparePDG(part->PdgCode(),signalMC->GetLegPDG(branch),signalMC->GetLegPDGexclude(branch),signalMC->GetCheckBothChargesLegs(branch))) return kFALSE;
1116 if(!CheckParticleSource(label, signalMC->GetLegSource(branch))) return kFALSE;
1119 AliVParticle* mcMother=0x0;
1121 if(signalMC->GetMotherPDG(branch)!=0 || signalMC->GetMotherSource(branch)!=AliDielectronSignalMC::kDontCare) {
1123 mLabel = GetMothersLabel(label);
1124 mcMother = GetMCTrackFromMCEvent(mLabel);
1126 if(!mcMother && !signalMC->GetMotherPDGexclude(branch)) return kFALSE;
1128 if(!ComparePDG((mcMother ? mcMother->PdgCode() : 0),signalMC->GetMotherPDG(branch),signalMC->GetMotherPDGexclude(branch),signalMC->GetCheckBothChargesMothers(branch))) return kFALSE;
1129 if(!CheckParticleSource(mLabel, signalMC->GetMotherSource(branch))) return kFALSE;
1131 //check for radiative deday
1132 if (!CheckRadiativeDecision(mLabel, signalMC)) return kFALSE;
1135 // check the grandmother
1136 AliVParticle* mcGrandMother=0x0;
1137 if(signalMC->GetGrandMotherPDG(branch)!=0 || signalMC->GetGrandMotherSource(branch)!=AliDielectronSignalMC::kDontCare) {
1140 gmLabel = GetMothersLabel(mLabel);
1141 mcGrandMother = static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(gmLabel));
1143 if(!mcGrandMother && !signalMC->GetGrandMotherPDGexclude(branch)) return kFALSE;
1145 if(!ComparePDG((mcGrandMother ? mcGrandMother->PdgCode() : 0),signalMC->GetGrandMotherPDG(branch),signalMC->GetGrandMotherPDGexclude(branch),signalMC->GetCheckBothChargesGrandMothers(branch))) return kFALSE;
1146 if(!CheckParticleSource(gmLabel, signalMC->GetGrandMotherSource(branch))) return kFALSE;
1154 //________________________________________________________________________________
1155 Bool_t AliDielectronMC::IsMCTruth(const AliDielectronPair* pair, const AliDielectronSignalMC* signalMC) const {
1157 // Check if the pair corresponds to the MC truth in signalMC
1161 const AliVParticle * mcD1 = pair->GetFirstDaughterP();
1162 const AliVParticle * mcD2 = pair->GetSecondDaughterP();
1163 Int_t labelD1 = (mcD1 ? TMath::Abs(mcD1->GetLabel()) : -1);
1164 Int_t labelD2 = (mcD2 ? TMath::Abs(mcD2->GetLabel()) : -1);
1165 Int_t d1Pdg = GetPdgFromLabel(labelD1);
1166 Int_t d2Pdg = GetPdgFromLabel(labelD2);
1169 AliVParticle* mcM1=0x0;
1170 AliVParticle* mcM2=0x0;
1173 AliVParticle* mcG1 = 0x0;
1174 AliVParticle* mcG2 = 0x0;
1176 // make direct(1-1 and 2-2) and cross(1-2 and 2-1) comparisons for the whole branch
1177 Bool_t directTerm = kTRUE;
1179 directTerm = directTerm && mcD1 && ComparePDG(d1Pdg,signalMC->GetLegPDG(1),signalMC->GetLegPDGexclude(1),signalMC->GetCheckBothChargesLegs(1))
1180 && CheckParticleSource(labelD1, signalMC->GetLegSource(1));
1182 directTerm = directTerm && mcD2 && ComparePDG(d2Pdg,signalMC->GetLegPDG(2),signalMC->GetLegPDGexclude(2),signalMC->GetCheckBothChargesLegs(2))
1183 && CheckParticleSource(labelD2, signalMC->GetLegSource(2));
1187 if(signalMC->GetMotherPDG(1)!=0 || signalMC->GetMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
1188 labelM1 = GetMothersLabel(labelD1);
1189 if(labelD1>-1 && labelM1>-1) mcM1 = GetMCTrackFromMCEvent(labelM1);
1190 directTerm = directTerm && (mcM1 || signalMC->GetMotherPDGexclude(1))
1191 && ComparePDG((mcM1 ? mcM1->PdgCode() : 0),signalMC->GetMotherPDG(1),signalMC->GetMotherPDGexclude(1),signalMC->GetCheckBothChargesMothers(1))
1192 && CheckParticleSource(labelM1, signalMC->GetMotherSource(1))
1193 && CheckRadiativeDecision(labelM1,signalMC);
1197 if(signalMC->GetMotherPDG(2)!=0 || signalMC->GetMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
1198 labelM2 = GetMothersLabel(labelD2);
1199 if(labelD2>-1 && labelM2>-1) mcM2 = GetMCTrackFromMCEvent(labelM2);
1200 directTerm = directTerm && (mcM2 || signalMC->GetMotherPDGexclude(2))
1201 && ComparePDG((mcM2 ? mcM2->PdgCode() : 0),signalMC->GetMotherPDG(2),signalMC->GetMotherPDGexclude(2),signalMC->GetCheckBothChargesMothers(2))
1202 && CheckParticleSource(labelM2, signalMC->GetMotherSource(2))
1203 && CheckRadiativeDecision(labelM2,signalMC);
1208 if(signalMC->GetGrandMotherPDG(1)!=0 || signalMC->GetGrandMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
1209 labelG1 = GetMothersLabel(labelM1);
1210 if(mcM1 && labelG1>-1) mcG1 = GetMCTrackFromMCEvent(labelG1);
1211 directTerm = directTerm && (mcG1 || signalMC->GetGrandMotherPDGexclude(1))
1212 && ComparePDG((mcG1 ? mcG1->PdgCode() : 0),signalMC->GetGrandMotherPDG(1),signalMC->GetGrandMotherPDGexclude(1),signalMC->GetCheckBothChargesGrandMothers(1))
1213 && CheckParticleSource(labelG1, signalMC->GetGrandMotherSource(1));
1217 if(signalMC->GetGrandMotherPDG(2)!=0 || signalMC->GetGrandMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
1218 labelG2 = GetMothersLabel(labelM2);
1219 if(mcM2 && labelG2>-1) mcG2 = GetMCTrackFromMCEvent(labelG2);
1220 directTerm = directTerm && (mcG2 || signalMC->GetGrandMotherPDGexclude(2))
1221 && ComparePDG((mcG2 ? mcG2->PdgCode() : 0),signalMC->GetGrandMotherPDG(2),signalMC->GetGrandMotherPDGexclude(2),signalMC->GetCheckBothChargesGrandMothers(2))
1222 && CheckParticleSource(labelG2, signalMC->GetGrandMotherSource(2));
1226 Bool_t crossTerm = kTRUE;
1228 crossTerm = crossTerm && mcD2
1229 && ComparePDG(d2Pdg,signalMC->GetLegPDG(1),signalMC->GetLegPDGexclude(1),signalMC->GetCheckBothChargesLegs(1))
1230 && CheckParticleSource(labelD2, signalMC->GetLegSource(1));
1232 crossTerm = crossTerm && mcD1
1233 && ComparePDG(d1Pdg,signalMC->GetLegPDG(2),signalMC->GetLegPDGexclude(2),signalMC->GetCheckBothChargesLegs(2))
1234 && CheckParticleSource(labelD1, signalMC->GetLegSource(2));
1237 if(signalMC->GetMotherPDG(1)!=0 || signalMC->GetMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
1238 if(!mcM2 && labelD2>-1) {
1239 labelM2 = GetMothersLabel(labelD2);
1240 if(labelM2>-1) mcM2 = GetMCTrackFromMCEvent(labelM2);
1242 crossTerm = crossTerm && (mcM2 || signalMC->GetMotherPDGexclude(1))
1243 && ComparePDG((mcM2 ? mcM2->PdgCode() : 0),signalMC->GetMotherPDG(1),signalMC->GetMotherPDGexclude(1),signalMC->GetCheckBothChargesMothers(1))
1244 && CheckParticleSource(labelM2, signalMC->GetMotherSource(1))
1245 && CheckRadiativeDecision(labelM2,signalMC);
1248 if(signalMC->GetMotherPDG(2)!=0 || signalMC->GetMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
1249 if(!mcM1 && labelD1>-1) {
1250 labelM1 = GetMothersLabel(labelD1);
1251 if(labelM1>-1) mcM1 = GetMCTrackFromMCEvent(labelM1);
1253 crossTerm = crossTerm && (mcM1 || signalMC->GetMotherPDGexclude(2))
1254 && ComparePDG((mcM1 ? mcM1->PdgCode() : 0),signalMC->GetMotherPDG(2),signalMC->GetMotherPDGexclude(2),signalMC->GetCheckBothChargesMothers(2))
1255 && CheckParticleSource(labelM1, signalMC->GetMotherSource(2))
1256 && CheckRadiativeDecision(labelM1,signalMC);
1260 if(signalMC->GetGrandMotherPDG(1)!=0 || signalMC->GetGrandMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
1262 labelG2 = GetMothersLabel(labelM2);
1263 if(labelG2>-1) mcG2 = GetMCTrackFromMCEvent(labelG2);
1265 crossTerm = crossTerm && (mcG2 || signalMC->GetGrandMotherPDGexclude(1))
1266 && ComparePDG((mcG2 ? mcG2->PdgCode() : 0),signalMC->GetGrandMotherPDG(1),signalMC->GetGrandMotherPDGexclude(1),signalMC->GetCheckBothChargesGrandMothers(1))
1267 && CheckParticleSource(labelG2, signalMC->GetGrandMotherSource(1));
1270 if(signalMC->GetGrandMotherPDG(2)!=0 || signalMC->GetGrandMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
1272 labelG1 = GetMothersLabel(labelM1);
1273 if(labelG1>-1) mcG1 = GetMCTrackFromMCEvent(labelG1);
1275 crossTerm = crossTerm && (mcG1 || signalMC->GetGrandMotherPDGexclude(2))
1276 && ComparePDG((mcG1 ? mcG1->PdgCode() : 0),signalMC->GetGrandMotherPDG(2),signalMC->GetGrandMotherPDGexclude(2),signalMC->GetCheckBothChargesGrandMothers(2))
1277 && CheckParticleSource(labelG1, signalMC->GetGrandMotherSource(2));
1280 Bool_t motherRelation = kTRUE;
1281 if(signalMC->GetMothersRelation()==AliDielectronSignalMC::kSame) {
1282 motherRelation = motherRelation && HaveSameMother(pair);
1284 if(signalMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent) {
1285 motherRelation = motherRelation && !HaveSameMother(pair);
1288 return ((directTerm || crossTerm) && motherRelation);
1293 //____________________________________________________________
1294 Bool_t AliDielectronMC::HaveSameMother(const AliDielectronPair * pair) const
1297 // Check whether two particles have the same mother
1300 const AliVParticle * daughter1 = pair->GetFirstDaughterP();
1301 const AliVParticle * daughter2 = pair->GetSecondDaughterP();
1302 if (!daughter1 || !daughter2) return 0;
1304 AliVParticle *mcDaughter1=GetMCTrackFromMCEvent(daughter1->GetLabel());
1305 AliVParticle *mcDaughter2=GetMCTrackFromMCEvent(daughter2->GetLabel());
1306 if (!mcDaughter1 || !mcDaughter2) return 0;
1308 Int_t labelMother1=-1;
1309 Int_t labelMother2=-1;
1311 if (mcDaughter1->IsA()==AliMCParticle::Class()){
1312 labelMother1=(static_cast<AliMCParticle*>(mcDaughter1))->GetMother();
1313 labelMother2=(static_cast<AliMCParticle*>(mcDaughter2))->GetMother();
1314 } else if (mcDaughter1->IsA()==AliAODMCParticle::Class()) {
1315 labelMother1=(static_cast<AliAODMCParticle*>(mcDaughter1))->GetMother();
1316 labelMother2=(static_cast<AliAODMCParticle*>(mcDaughter2))->GetMother();
1319 Bool_t sameMother=(labelMother1>-1)&&(labelMother2>-1)&&(labelMother1==labelMother2);
1324 //________________________________________________________________
1325 Int_t AliDielectronMC::IsJpsiPrimary(const AliDielectronPair * pair)
1327 // return: "0" for primary jpsi
1328 // "1" for secondary jpsi (from beauty)
1329 // "2" for background
1330 if(!HaveSameMother(pair)) return 2;
1331 AliVParticle *mcDaughter1=GetMCTrackFromMCEvent((pair->GetFirstDaughterP())->GetLabel());
1332 Int_t labelMother=-1;
1334 if (mcDaughter1->IsA()==AliMCParticle::Class()){
1335 labelMother=(static_cast<AliMCParticle*>(mcDaughter1))->GetMother();
1336 } else if (mcDaughter1->IsA()==AliAODMCParticle::Class()) {
1337 labelMother=(static_cast<AliAODMCParticle*>(mcDaughter1))->GetMother();
1340 AliVParticle* mcMother=GetMCTrackFromMCEvent(labelMother);
1341 if(!IsMCMotherToEE(mcMother,443)) return 2;
1342 return IsJpsiPrimary(mcMother);
1345 //______________________________________________________________
1346 Int_t AliDielectronMC::IsJpsiPrimary(const AliVParticle * particle)
1348 // return: "0" for primary jpsi
1349 // "1" for secondary jpsi (come from B decay)
1353 if (particle->IsA()==AliMCParticle::Class()){
1354 labelMoth = (static_cast<const AliMCParticle*>(particle))->GetMother();
1356 particle = GetMCTrackFromMCEvent(labelMoth);
1357 pdgCode = TMath::Abs((static_cast<const AliMCParticle*>(particle))->PdgCode());
1358 if((pdgCode>500 && pdgCode<600) || (pdgCode>5000 && pdgCode<6000)) return 1;
1359 labelMoth = (static_cast<const AliMCParticle*>(particle))->GetMother();
1362 else if (particle->IsA()==AliAODMCParticle::Class()){
1363 labelMoth = (static_cast<const AliAODMCParticle*>(particle))->GetMother();
1365 particle = GetMCTrackFromMCEvent(labelMoth);
1366 pdgCode = TMath::Abs((static_cast<const AliAODMCParticle*>(particle))->PdgCode());
1367 if((pdgCode>500 && pdgCode<600) || (pdgCode>5000 && pdgCode<6000)) return 1;
1368 labelMoth = (static_cast<const AliAODMCParticle*>(particle))->GetMother();
1375 Bool_t AliDielectronMC::GetPrimaryVertex(Double_t &primVtxX, Double_t &primVtxY, Double_t &primVtxZ){
1377 if(fAnaType == kESD){
1378 const AliVVertex* mcVtx = fMCEvent->GetPrimaryVertex();
1379 if(!mcVtx) return kFALSE;
1380 primVtxX = mcVtx->GetX();
1381 primVtxY = mcVtx->GetY();
1382 primVtxZ = mcVtx->GetZ();
1383 }else if(fAnaType == kAOD){
1384 AliAODEvent *aod=((AliAODInputHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()))->GetEvent();
1385 if(!aod) return kFALSE;
1386 AliAODMCHeader *mcHead = dynamic_cast<AliAODMCHeader*>(aod->FindListObject(AliAODMCHeader::StdBranchName()));
1387 if(!mcHead) return kFALSE;
1388 primVtxX = mcHead->GetVtxX();
1389 primVtxY = mcHead->GetVtxY();
1390 primVtxZ = mcHead->GetVtxZ();