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