b4e37bda51f721c1ef207c637f8d1f962820457d
[u/mrichter/AliRoot.git] / PWG3 / 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 /* $Id$ */
17
18 //#####################################################
19 //#                                                   # 
20 //#              Class AliDielectronMC                #
21 //#       Cut Class for Jpsi->e+e- analysis           #
22 //#                                                   #
23 //#   by WooJin J. Park, GSI / W.J.Park@gsi.de        #
24 //#                                                   #
25 //#####################################################
26
27 #include <AliAnalysisManager.h>
28 #include <AliAODHandler.h>
29 #include <AliESDInputHandler.h>
30 #include <AliAODInputHandler.h>
31 #include <AliMCEventHandler.h>
32 #include <AliMCEvent.h>
33 #include <AliMCParticle.h>
34 #include <AliAODMCParticle.h>
35 #include <AliStack.h>
36 #include <AliESDEvent.h>
37 #include <AliAODEvent.h>
38 #include <AliESDtrack.h>
39 #include <AliAODTrack.h>
40 #include <AliLog.h>
41
42 #include <TClonesArray.h>
43 #include "AliDielectronMC.h"
44
45 AliDielectronMC* AliDielectronMC::fgInstance=0x0;
46
47 //____________________________________________________________
48 AliDielectronMC* AliDielectronMC::Instance()
49 {
50   //
51   // return pointer to singleton implementation
52   //
53   if (fgInstance) return fgInstance;
54
55   AnalysisType type=kUNSET;
56   Bool_t hasMC=kFALSE;
57   if (AliAnalysisManager::GetAnalysisManager()){
58     if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliESDInputHandler::Class()) type=kESD;
59     else if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliAODInputHandler::Class()) type=kAOD;
60
61     AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
62     if(type == kESD) hasMC=mcHandler!=0x0;
63     }
64  
65   fgInstance=new AliDielectronMC(type);
66  
67   fgInstance->SetHasMC(hasMC);
68    
69   return fgInstance;
70 }
71
72 //____________________________________________________________
73 AliDielectronMC::AliDielectronMC(AnalysisType type):
74   fMCEvent(0x0),
75   fStack(0x0),
76   fAnaType(type),
77   fHasMC(kTRUE),
78   fMcArray(0x0)
79 {
80   //
81   // default constructor
82   //
83 }
84
85
86 //____________________________________________________________
87 AliDielectronMC::~AliDielectronMC()
88 {
89   //
90   // default destructor
91   //
92   
93 }
94
95 //____________________________________________________________
96 void AliDielectronMC::Initialize()
97 {
98   //
99   // initialize MC class
100   //
101   if (!ConnectMCEvent()) AliError("Initialization of MC object failed!");
102 }
103
104 //____________________________________________________________
105 Int_t AliDielectronMC::GetNMCTracks()
106 {
107   //
108   //  return the number of generated tracks from MC event
109   //
110   if(fAnaType == kESD){
111     if (!fMCEvent){ AliError("No fMCEvent"); return 0; }
112     return fMCEvent->GetNumberOfTracks();}
113   else if(fAnaType == kAOD){
114     if(!fMcArray) { AliError("No fMcArray"); return 0; }
115     return fMcArray->GetEntriesFast();
116   }
117   return 0;
118 }
119
120 //____________________________________________________________
121 Int_t AliDielectronMC::GetNMCTracksFromStack()
122 {
123   //
124   //  return the number of generated tracks from stack
125   //
126   if (!fStack){ AliError("No fStack"); return -999; }
127   return fStack->GetNtrack();
128 }
129
130 //____________________________________________________________
131 Int_t AliDielectronMC::GetNPrimary()
132 {
133   //
134   //  return the number of primary track from MC event
135   //
136   if (!fMCEvent){ AliError("No fMCEvent"); return 0; }
137   return fMCEvent->GetNumberOfPrimaries();
138 }
139
140 //____________________________________________________________
141 Int_t AliDielectronMC::GetNPrimaryFromStack()
142 {
143   //
144   //  return the number of primary track from stack
145   //
146   if (!fStack){ AliError("No fStack"); return -999; }
147   return fStack->GetNprimary();
148 }
149
150 //____________________________________________________________
151 AliVParticle* AliDielectronMC::GetMCTrackFromMCEvent(Int_t itrk)
152 {
153   //
154   // return MC track directly from MC event
155   //
156   if (itrk<0) return NULL;
157   AliVParticle * track=0x0;
158   if(fAnaType == kESD){
159     if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
160     track = fMCEvent->GetTrack(itrk); //  tracks from MC event (ESD)
161   } else if(fAnaType == kAOD) {
162     if (!fMcArray){ AliError("No fMCEvent"); return NULL;}
163     track = (AliVParticle*)fMcArray->At(itrk); //  tracks from MC event (AOD)
164   }
165   return track;
166 }
167
168 //____________________________________________________________
169 Bool_t AliDielectronMC::ConnectMCEvent()
170 {
171   //
172   // connect stack object from the mc handler
173   //
174   if(fAnaType == kESD){
175     fMCEvent=0x0;
176     AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
177     if (!mcHandler){ AliError("Could not retrive MC event handler!"); return kFALSE; }
178     if (!mcHandler->InitOk() ) return kFALSE;
179     if (!mcHandler->TreeK() )  return kFALSE;
180     if (!mcHandler->TreeTR() ) return kFALSE;
181     
182     AliMCEvent* mcEvent = mcHandler->MCEvent();
183     if (!mcEvent){ AliError("Could not retrieve MC event!"); return kFALSE; }
184     fMCEvent = mcEvent;
185     
186     if (!UpdateStack()) return kFALSE;
187   }
188   else if(fAnaType == kAOD)
189   {
190     fMcArray = 0x0;
191     AliAODEvent *aod=((AliAODInputHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()))->GetEvent();
192     fMcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
193     if(!fMcArray) return kFALSE;
194   }
195   return kTRUE;
196 }
197
198 //____________________________________________________________
199 Bool_t AliDielectronMC::UpdateStack()
200 {
201   //
202   // update stack with new event
203   //
204   if (!fMCEvent){ AliError("No fMCEvent"); return kFALSE;}
205   AliStack* stack = fMCEvent->Stack();
206   if (!stack){ AliError("Could not retrive stack!"); return kFALSE; }
207   fStack = stack;
208   return kTRUE;
209 }
210
211 //____________________________________________________________
212 AliMCParticle* AliDielectronMC::GetMCTrack( const AliESDtrack* _track)
213 {
214   //
215   // return MC track
216   //
217   if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
218
219   Int_t nStack = fMCEvent->GetNumberOfTracks();
220   Int_t label = TMath::Abs(_track->GetLabel());
221   if(label>nStack)return NULL;
222
223   AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(label));
224   return mctrack;
225 }
226
227
228 //____________________________________________________________
229 AliAODMCParticle* AliDielectronMC::GetMCTrack( const AliAODTrack* _track)
230 {
231   //
232   // return MC track
233   //
234  if(!fMcArray) { AliError("No fMCArray"); return NULL;}
235  Int_t label = _track->GetLabel();
236  if(label < 0) return NULL;
237  AliAODMCParticle *mctrack = (AliAODMCParticle*)fMcArray->At(label);
238  return mctrack; 
239 }
240
241 //____________________________________________________________
242 TParticle* AliDielectronMC::GetMCTrackFromStack(const AliESDtrack* _track)
243 {
244   //
245   // return MC track from stack
246   //
247   Int_t label = TMath::Abs(_track->GetLabel());
248   if (!fStack) AliWarning("fStack is not available. Update stack first.");
249   TParticle* mcpart = fStack->Particle(label);
250   if (!mcpart) return NULL;
251   return mcpart;
252 }
253
254 //____________________________________________________________
255 AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliESDtrack* _track)
256 {
257   //
258   // return MC track mother
259   //
260   AliMCParticle* mcpart = GetMCTrack(_track);
261   if (!mcpart) return NULL;
262   AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(mcpart->GetMother()));
263   if (!mcmother) return NULL;
264   return mcmother;
265 }
266
267 //______________________________________________________________
268 AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODTrack* _track)
269 {
270  //
271  // return MC track mother
272  //
273  AliAODMCParticle* mcpart = GetMCTrack(_track);
274  if (!mcpart) return NULL;
275  if(mcpart->GetMother() < 0) return NULL;
276  AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMcArray->At(mcpart->GetMother()));
277  if (!mcmother) return NULL;
278  return mcmother;
279 }
280 //____________________________________________________________
281 AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliMCParticle* _particle){
282   //
283   // return MC track mother
284   //
285   AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(_particle->GetMother()));
286   return mcmother;
287 }
288
289 //____________________________________________________________
290 AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODMCParticle* _particle){
291   //
292   // return MC track mother
293   //
294   if( _particle->GetMother() < 0) return NULL;
295   AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMcArray->At(_particle->GetMother()));
296   return mcmother;
297 }
298
299 //____________________________________________________________
300 TParticle* AliDielectronMC::GetMCTrackMotherFromStack(const AliESDtrack* _track)
301 {
302   //
303   // return MC track mother from stack
304   //
305   TParticle* mcpart = GetMCTrackFromStack(_track);
306   if ( !mcpart || mcpart->GetFirstMother()<=0 ) return NULL; 
307   TParticle* mcmother = fStack->Particle(mcpart->GetFirstMother());
308   if (!mcmother) return NULL;
309   return mcmother;
310 }
311
312 //____________________________________________________________
313 Int_t AliDielectronMC::GetMCPID(const AliESDtrack* _track)
314 {
315   //
316   // return PDG code of the track from the MC truth info
317   //
318   AliMCParticle* mcpart = GetMCTrack(_track);
319   if (!mcpart) return -999;
320   return mcpart->PdgCode();
321 }
322
323 //__________________________________________________________
324 Int_t AliDielectronMC::GetMCPID(const AliAODTrack* _track)
325 {
326  //
327  // return PDG code of the track from the MC truth info
328  //
329  AliAODMCParticle* mcpart = GetMCTrack(_track);
330  if (!mcpart) return -999;
331  return mcpart->PdgCode();
332 }
333
334 //____________________________________________________________
335 Int_t AliDielectronMC::GetMCPIDFromStack(const AliESDtrack* _track)
336 {
337   // 
338   // return MC PDG code from stack
339   //
340   TParticle* mcpart = GetMCTrackFromStack(_track);
341   if (!mcpart) return -999;
342   return mcpart->GetPdgCode();
343 }
344
345 //____________________________________________________________
346 Int_t AliDielectronMC::GetMotherPDG( const AliESDtrack* _track)
347 {
348   //
349   // return PDG code of the mother track from the MC truth info
350   //
351   AliMCParticle* mcmother = GetMCTrackMother(_track);
352   if (!mcmother) return -999;
353   return mcmother->PdgCode();
354 }
355
356 //________________________________________________________
357 Int_t AliDielectronMC::GetMotherPDG( const AliAODTrack* _track)
358 {
359   //
360   // return PDG code of the mother track from the MC truth info
361   //
362   AliAODMCParticle* mcmother = GetMCTrackMother(_track);
363   if (!mcmother) return -999;
364   return mcmother->PdgCode();
365 }
366
367 //____________________________________________________________
368 Int_t AliDielectronMC::GetMotherPDGFromStack(const AliESDtrack* _track)
369 {
370   //
371   // return PDG code of the mother track from stack
372   //
373   TParticle* mcmother = GetMCTrackMotherFromStack(_track);
374   if (!mcmother) return -999;
375   return mcmother->GetPdgCode();
376 }
377
378 //____________________________________________________________
379 Int_t AliDielectronMC::GetMCProcess(const AliESDtrack* _track)
380 {
381   //
382   // return process number of the track
383   //
384   AliMCParticle* mcpart = GetMCTrack(_track);
385   if (!mcpart) return -999;
386   return 0;
387 }
388
389 //____________________________________________________________
390 Int_t AliDielectronMC::GetMCProcessFromStack(const AliESDtrack* _track)
391 {
392   //
393   // return process number of the track
394   //
395   TParticle* mcpart = GetMCTrackFromStack(_track);
396   if (!mcpart) return -999;
397   return mcpart->GetUniqueID();
398 }
399
400 //____________________________________________________________
401 Int_t AliDielectronMC::NumberOfDaughters(const AliESDtrack* track)
402 {
403   //
404   // returns the number of daughters
405   //
406   AliMCParticle *mcmother=GetMCTrackMother(track);
407   if(!mcmother||!mcmother->Particle()) return -999;
408 //   return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
409   return mcmother->Particle()->GetNDaughters();
410 }
411
412 //_________________________________________________________
413 Int_t AliDielectronMC::NumberOfDaughters(const AliAODTrack* track)
414 {
415   //
416   // returns the number of daughters
417   //
418   AliAODMCParticle *mcmother=GetMCTrackMother(track);
419   if(!mcmother) return -999;
420   return NumberOfDaughters(mcmother);  
421
422 }
423
424 //____________________________________________________________
425 Int_t AliDielectronMC::NumberOfDaughters(const AliMCParticle* particle)
426 {
427   //
428   // returns the number of daughters
429   //
430   AliMCParticle *mcmother=GetMCTrackMother(particle);
431   if(!mcmother||!mcmother->Particle()) return -999;
432   //return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
433   return mcmother->Particle()->GetNDaughters();
434 }
435
436 //____________________________________________________________
437 Int_t AliDielectronMC::NumberOfDaughters(const AliAODMCParticle* particle)
438 {
439   //
440   // returns the number of daughters
441   //
442   AliAODMCParticle *mcmother=GetMCTrackMother(particle);
443   if(!mcmother) return -999;
444   return mcmother->GetNDaughters();
445 }
446
447 //____________________________________________________________
448 Int_t AliDielectronMC::GetMCProcessMother(const AliESDtrack* _track)
449 {
450   //
451   // return process number of the mother of the track
452   //
453   AliMCParticle* mcmother = GetMCTrackMother(_track);
454   if (!mcmother) return -999;
455   return 0;
456 }
457
458 //____________________________________________________________
459 Int_t AliDielectronMC::GetMCProcessMotherFromStack(const AliESDtrack* _track)
460 {
461   //
462   // return process number of the mother of the track
463   //
464   TParticle* mcmother = GetMCTrackMotherFromStack(_track);
465   if (!mcmother) return -999;
466   return mcmother->GetUniqueID();
467 }
468
469 //____________________________________________________________
470 Bool_t AliDielectronMC::IsMCMotherToEE(const AliVParticle *particle, Int_t pdgMother)
471 {
472   //
473   // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
474   //
475   if (fAnaType==kESD && !fMCEvent) return kFALSE;
476   if (fAnaType==kAOD && !fMcArray) return kFALSE;
477   if (!particle) return kFALSE;
478   
479   if (particle->IsA()==AliMCParticle::Class()){
480     return IsMCMotherToEEesd(static_cast<const AliMCParticle*>(particle),pdgMother);
481   } else if (particle->IsA()==AliAODMCParticle::Class()){
482    return IsMCMotherToEEaod(static_cast<const AliAODMCParticle*>(particle),pdgMother);
483   } else {
484     AliError("Unknown particle type");
485   }
486   return kFALSE;
487 }
488
489 //____________________________________________________________
490 Bool_t AliDielectronMC::IsMCMotherToEEesd(const AliMCParticle *particle, Int_t pdgMother)
491 {
492   //
493   // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
494   // ESD case
495   //
496   
497   //check pdg code
498   if (particle->PdgCode()!=pdgMother) return kFALSE;
499   Int_t ifirst = particle->GetFirstDaughter();
500   Int_t ilast  = particle->GetLastDaughter();
501   
502   //check number of daughters
503   if ((ilast-ifirst)!=1) return kFALSE;
504   AliMCParticle *firstD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ifirst));
505   AliMCParticle *secondD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ilast));
506
507   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
508   if (firstD->Charge()>0){
509     if (firstD->PdgCode()!=-11) return kFALSE;
510     if (secondD->PdgCode()!=11) return kFALSE;
511   }else{
512     if (firstD->PdgCode()!=11) return kFALSE;
513     if (secondD->PdgCode()!=-11) return kFALSE;
514   }
515   
516   return kTRUE;
517 }
518
519 //____________________________________________________________
520 Bool_t AliDielectronMC::IsMCMotherToEEaod(const AliAODMCParticle *particle, Int_t pdgMother)
521 {
522   //
523   // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
524   // AOD case
525   //
526
527   if (particle->GetPdgCode()!=pdgMother) return kFALSE;
528   if (particle->GetNDaughters()!=2) return kFALSE;
529  
530   Int_t ifirst = particle->GetDaughter(0);
531   Int_t ilast  = particle->GetDaughter(1);
532   
533   //check number of daughters
534   if ((ilast-ifirst)!=1) return kFALSE;
535   
536   AliAODMCParticle *firstD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ifirst));
537   AliAODMCParticle *secondD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ilast));
538    
539   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
540  
541   if (firstD->Charge()>0){
542     if (firstD->GetPdgCode()!=-11) return kFALSE;
543     if (secondD->GetPdgCode()!=11) return kFALSE;
544   }else{
545     if (firstD->GetPdgCode()!=11) return kFALSE;
546     if (secondD->GetPdgCode()!=-11) return kFALSE;
547   }
548   return kTRUE;
549 }
550
551 //____________________________________________________________
552 Int_t AliDielectronMC::GetLabelMotherWithPdg(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
553 {
554   //
555   // test if mother of particle 1 and 2 has pdgCode pdgMother and is the same;
556   //
557   if (fAnaType==kESD){
558   if (!fMCEvent) return -1;
559   return GetLabelMotherWithPdgESD(particle1, particle2, pdgMother);
560   }
561   else if (fAnaType==kAOD)
562   {
563   if (!fMcArray) return -1;
564   return GetLabelMotherWithPdgAOD(particle1, particle2, pdgMother);
565   }  
566
567   return -1;
568 }
569
570 //____________________________________________________________
571 Int_t AliDielectronMC::GetLabelMotherWithPdgESD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
572 {
573   //
574   // test if mother of particle 1 and 2 has pdgCode +-11 (electron),
575   //    have the same mother and the mother had pdg code pdgMother
576   // ESD case
577   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
578   //
579   
580   AliMCParticle *mcPart1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(particle1->GetLabel()));
581   AliMCParticle *mcPart2=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(particle2->GetLabel()));
582   
583   if (!mcPart1||!mcPart2) return -1;
584   
585   Int_t lblMother1=mcPart1->GetMother();
586   Int_t lblMother2=mcPart2->GetMother();
587   
588   AliMCParticle *mcMother1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(lblMother1));
589   if (!mcMother1) return -1;
590   if (lblMother1!=lblMother2) return -1;
591   if (TMath::Abs(mcPart1->PdgCode())!=11) return -1;
592   if (mcPart1->PdgCode()!=-mcPart2->PdgCode()) return -1;
593   if (mcMother1->PdgCode()!=pdgMother) return -1;
594   
595   return lblMother1;
596 }
597
598 //____________________________________________________________
599 Int_t AliDielectronMC::GetLabelMotherWithPdgAOD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
600 {
601   //
602   // test if mother of particle 1 and 2 has pdgCode +-11 (electron),
603   //    have the same mother and the mother had pdg code pdgMother
604   // AOD case
605   //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
606   //
607   AliAODMCParticle *mcPart1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(particle1->GetLabel()));
608   AliAODMCParticle *mcPart2=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(particle2->GetLabel()));
609   
610   if (!mcPart1||!mcPart2) return -1;
611   
612   Int_t lblMother1=mcPart1->GetMother();
613   Int_t lblMother2=mcPart2->GetMother();
614   
615   AliAODMCParticle *mcMother1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(lblMother1));
616   
617   if (!mcMother1) return -1;
618   if (lblMother1!=lblMother2) return -1;
619   if (TMath::Abs(mcPart1->GetPdgCode())!=11) return -1;
620   if (mcPart1->GetPdgCode()!=-mcPart2->GetPdgCode()) return -1;
621   if (mcMother1->GetPdgCode()!=pdgMother) return -1;
622   
623   return lblMother1;
624 }
625
626 //____________________________________________________________
627 void AliDielectronMC::GetDaughters(const TObject *mother, AliVParticle* &d1, AliVParticle* &d2)
628 {
629   //
630   // Get First two daughters of the mother
631   //
632   Int_t lblD1=-1;
633   Int_t lblD2=-1;
634   d1=0;
635   d2=0;
636   if (fAnaType==kAOD){
637     if(!fMcArray) return;
638     const AliAODMCParticle *aodMother=static_cast<const AliAODMCParticle*>(mother);
639     lblD1=aodMother->GetDaughter(0);
640     lblD2=aodMother->GetDaughter(1);
641     d1 = (AliVParticle*)fMcArray->At(lblD1);
642     d2 = (AliVParticle*)fMcArray->At(lblD2);
643    } else if (fAnaType==kESD){
644     if (!fMCEvent) return;
645     const AliMCParticle *aodMother=static_cast<const AliMCParticle*>(mother);
646     lblD1=aodMother->GetFirstDaughter();
647     lblD2=aodMother->GetLastDaughter();
648     d1=fMCEvent->GetTrack(lblD1);
649     d2=fMCEvent->GetTrack(lblD2);
650    }
651 }
652
653 //____________________________________________________________
654 Bool_t AliDielectronMC::HaveSameMother(const AliDielectronPair * pair)
655 {
656   //
657   // Check whether two particles have the same mother
658   //
659
660   const AliVParticle * daughter1 = pair->GetFirstDaughter();
661   const AliVParticle * daughter2 = pair->GetSecondDaughter();
662
663   AliVParticle *mcDaughter1=GetMCTrackFromMCEvent(daughter1->GetLabel());
664   AliVParticle *mcDaughter2=GetMCTrackFromMCEvent(daughter2->GetLabel());
665   if (!mcDaughter1 || !mcDaughter2) return 0;
666
667   Int_t labelMother1=-1;
668   Int_t labelMother2=-1;
669
670   if (mcDaughter1->IsA()==AliMCParticle::Class()){
671     labelMother1=(static_cast<AliMCParticle*>(mcDaughter1))->GetMother();
672     labelMother2=(static_cast<AliMCParticle*>(mcDaughter2))->GetMother();
673   } else if (mcDaughter1->IsA()==AliAODMCParticle::Class()) {
674     labelMother1=(static_cast<AliAODMCParticle*>(mcDaughter1))->GetMother();
675     labelMother2=(static_cast<AliAODMCParticle*>(mcDaughter2))->GetMother();
676   }
677
678   Bool_t sameMother=(labelMother1>-1)&&(labelMother2>-1)&&(labelMother1==labelMother2);
679
680   return sameMother;
681 }
682
683 //________________________________________________________________
684 Int_t AliDielectronMC::IsJpsiPrimary(const AliDielectronPair * pair)
685 {
686  // return: "0" for primary jpsi 
687  //         "1" for secondary jpsi (from beauty)
688  //         "2" for background  
689  if(!HaveSameMother(pair)) return 2;
690  AliVParticle *mcDaughter1=GetMCTrackFromMCEvent((pair->GetFirstDaughter())->GetLabel());
691  Int_t labelMother=-1;
692
693   if (mcDaughter1->IsA()==AliMCParticle::Class()){
694      labelMother=(static_cast<AliMCParticle*>(mcDaughter1))->GetMother();
695      } else if (mcDaughter1->IsA()==AliAODMCParticle::Class()) {
696      labelMother=(static_cast<AliAODMCParticle*>(mcDaughter1))->GetMother();
697      }
698
699  AliVParticle* mcMother=GetMCTrackFromMCEvent(labelMother);
700  if(!IsMCMotherToEE(mcMother,443)) return 2;
701  return IsJpsiPrimary(mcMother);
702 }
703
704 //______________________________________________________________
705 Int_t AliDielectronMC::IsJpsiPrimary(const AliVParticle * particle)
706 {
707   // return: "0" for primary jpsi
708   //         "1" for secondary jpsi (come from B decay)
709  Int_t labelMoth=-1;
710  Int_t pdgCode;
711
712  if (particle->IsA()==AliMCParticle::Class()){
713      labelMoth = (static_cast<const AliMCParticle*>(particle))->GetMother();
714      while(labelMoth>0){
715        particle = GetMCTrackFromMCEvent(labelMoth);
716        pdgCode = TMath::Abs((static_cast<const AliMCParticle*>(particle))->PdgCode());
717        if((pdgCode>500 && pdgCode<600) || (pdgCode>5000 && pdgCode<6000)) return 1;
718        labelMoth = (static_cast<const AliMCParticle*>(particle))->GetMother();
719        }
720     }
721  else if (particle->IsA()==AliAODMCParticle::Class()){
722      labelMoth = (static_cast<const AliAODMCParticle*>(particle))->GetMother();
723      while(labelMoth>0){
724      particle = GetMCTrackFromMCEvent(labelMoth);
725      pdgCode = TMath::Abs((static_cast<const AliAODMCParticle*>(particle))->PdgCode());
726      if((pdgCode>500 && pdgCode<600) || (pdgCode>5000 && pdgCode<6000)) return 1;
727      labelMoth = (static_cast<const AliAODMCParticle*>(particle))->GetMother();
728      }
729   }
730   return 0;
731 }