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