AliPhysics  59e0e03 (59e0e03)
AliAnalysisTaskSEDmesonsFilterCJ.cxx
Go to the documentation of this file.
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 // Task for selecting D mesons to be used as an input for D-Jet correlations
18 //
19 //-----------------------------------------------------------------------
20 // Authors:
21 // C. Bianchin (Utrecht University) chiara.bianchin@cern.ch
22 // A. Grelli (Utrecht University) a.grelli@uu.nl
23 // X. Zhang (LBNL) XMZhang@lbl.gov
24 // S. Aiola (Yale University) salvatore.aiola@cern.ch
25 // S. Antônio (University of São Paulo / Utrecht University) antonio.silva@cern.ch
26 // B. Trzeciak (Utrecht University) barbara.antonia.trzeciak@cern.ch
27 //-----------------------------------------------------------------------
28 
29 #include <TDatabasePDG.h>
30 #include <TParticle.h>
31 #include <TVector3.h>
32 #include <TROOT.h>
33 #include <TH3F.h>
34 #include <TRandom3.h>
35 
36 #include "AliLog.h"
38 #include "AliRDHFCutsD0toKpi.h"
39 #include "AliAODMCHeader.h"
40 #include "AliAODHandler.h"
41 #include "AliAnalysisManager.h"
42 #include "AliLog.h"
43 #include "AliAODVertex.h"
44 #include "AliAODRecoCascadeHF.h"
46 #include "AliESDtrack.h"
47 #include "AliAODMCParticle.h"
48 #include "AliEmcalParticle.h"
49 #include "AliParticleContainer.h"
50 #include "AliAnalysisDataSlot.h"
51 #include "AliAnalysisDataContainer.h"
52 #include "AliStack.h"
53 #include "AliAnalysisVertexingHF.h"
54 #include "AliVertexingHFUtils.h"
55 
57 
59 
60 //_______________________________________________________________________________
63  fUseMCInfo(kFALSE),
64  fBuildRMEff(kFALSE),
65  fUsePythia(kFALSE),
66  fUseRejTracks(kFALSE),
67  fTrackIneff(0),
68  fUseReco(kTRUE),
69  fCandidateType(0),
70  fCandidateName(""),
71  fPDGmother(0),
72  fNProngs(0),
73  fBranchName(""),
74  fCuts(0),
75  fMinMass(0.),
76  fMaxMass(0.),
77  fInhibitTask(kFALSE),
78  fCombineDmesons(kFALSE),
79  fMultCand(kFALSE),
80  fAnalyseCand(0),
81  fRejectQuarkNotFound(kTRUE),
82  fRejectDfromB(kTRUE),
83  fKeepOnlyDfromB(kFALSE),
84  fAodEvent(0),
85  fMCHeader(0),
86  fRan(0),
87  fArrayDStartoD0pi(0),
88  fMCarray(0),
89  fCandidateArray(0),
90  fSideBandArray(0),
91  fCombinedDmesons(0),
92  fCombinedDmesonsBkg(0),
93  fMCCombinedDmesons(0),
94  fNCand(0),
95  fNSBCand(0),
96  fHistStat(0),
97  fHistNSBCandEv(0),
98  fHistNCandEv(0),
99  fHistImpParS(0),
100  fHistImpParB(0),
101  fHistPtPion(0),
102  fHistInvMassPtD(0),
103  fHistInvMassS(0),
104  fHistInvMassB(0),
105  fHistAlphaDDS(0),
106  fHistAlphaDpisS(0),
107  fHistAlphaDpiS(0),
108  fHistAlphaDKS(0),
109  fHistAlphaDDB(0),
110  fHistAlphaDpisB(0),
111  fHistAlphaDpiB(0),
112  fHistAlphaDKB(0),
113  fHistDeltaRDDS(0),
114  fHistDeltaRDpisS(0),
115  fHistDeltaRDpiS(0),
116  fHistDeltaRDKS(0),
117  fHistDeltaRDDB(0),
118  fHistDeltaRDpisB(0),
119  fHistDeltaRDpiB(0),
120  fHistDeltaRDKB(0),
121  fHistAlphaDpiR(0),
122  fHistAlphaDKR(0),
123  fHistDeltaRDpiR(0),
124  fHistDeltaRDKR(0)
125 {
126  //
127  // Default constructor
128  //
129 
130  for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
131  for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
132 
133  fNeedEmcalGeom = kFALSE;
134 }
135 
136 //_______________________________________________________________________________
138  AliAnalysisTaskEmcal(name, kTRUE),
139  fUseMCInfo(kFALSE),
140  fBuildRMEff(kFALSE),
141  fUsePythia(kFALSE),
142  fUseRejTracks(kFALSE),
143  fTrackIneff(0),
144  fUseReco(kTRUE),
145  fCandidateType(candtype),
146  fCandidateName(""),
147  fPDGmother(0),
148  fNProngs(0),
149  fBranchName(""),
150  fCuts(cuts),
151  fMinMass(0.),
152  fMaxMass(0.),
153  fInhibitTask(kFALSE),
154  fCombineDmesons(kFALSE),
155  fMultCand(kFALSE),
156  fAnalyseCand(0),
157  fRejectQuarkNotFound(kTRUE),
158  fRejectDfromB(kTRUE),
159  fKeepOnlyDfromB(kFALSE),
160  fAodEvent(0),
161  fMCHeader(0),
162  fRan(0),
163  fArrayDStartoD0pi(0),
164  fMCarray(0),
165  fCandidateArray(0),
166  fSideBandArray(0),
167  fCombinedDmesons(0),
168  fCombinedDmesonsBkg(0),
169  fMCCombinedDmesons(0),
170  fNCand(0),
171  fNSBCand(0),
172  fHistStat(0),
173  fHistNSBCandEv(0),
174  fHistNCandEv(0),
175  fHistImpParS(0),
176  fHistImpParB(0),
177  fHistPtPion(0),
178  fHistInvMassPtD(0),
179  fHistInvMassS(0),
180  fHistInvMassB(0),
181  fHistAlphaDDS(0),
182  fHistAlphaDpisS(0),
183  fHistAlphaDpiS(0),
184  fHistAlphaDKS(0),
185  fHistAlphaDDB(0),
186  fHistAlphaDpisB(0),
187  fHistAlphaDpiB(0),
188  fHistAlphaDKB(0),
189  fHistDeltaRDDS(0),
190  fHistDeltaRDpisS(0),
191  fHistDeltaRDpiS(0),
192  fHistDeltaRDKS(0),
193  fHistDeltaRDDB(0),
194  fHistDeltaRDpisB(0),
195  fHistDeltaRDpiB(0),
196  fHistDeltaRDKB(0),
197  fHistAlphaDpiR(0),
198  fHistAlphaDKR(0),
199  fHistDeltaRDpiR(0),
200  fHistDeltaRDKR(0)
201 {
202  //
203  // Constructor. Initialization of Inputs and Outputs
204  //
205 
206  Info("AliAnalysisTaskSEDmesonsFilterCJ","Calling Constructor");
207 
208  fNeedEmcalGeom = kFALSE;
209 
210  for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
211  for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
212 
213  const Int_t nptbins = fCuts->GetNPtBins();
214  Float_t defaultSigmaD013[20]={0.012, 0.012, 0.012, 0.015, 0.015,0.018,0.018,0.020,0.020,0.030,0.030,0.037,0.040,0.040,0.040,0.040,0.040,0.040,0.040,0.040};
215 
216  switch (fCandidateType) {
217  case 0 :
218  fCandidateName = "D0";
219  fPDGmother = 421;
220  fNProngs = 2;
221  fPDGdaughters[0] = 211; // pi
222  fPDGdaughters[1] = 321; // K
223  fPDGdaughters[2] = 0; // empty
224  fPDGdaughters[3] = 0; // empty
225  fBranchName = "D0toKpi";
226  break;
227  case 1 :
228  fCandidateName = "DStar";
229  fPDGmother = 413;
230  fNProngs = 3;
231  fPDGdaughters[1] = 211; // pi soft
232  fPDGdaughters[0] = 421; // D0
233  fPDGdaughters[2] = 211; // pi fromD0
234  fPDGdaughters[3] = 321; // K from D0
235  fBranchName = "Dstar";
236 
237  if (nptbins<20) {
238  for (Int_t ipt=0; ipt<nptbins;ipt++) fSigmaD0[ipt] = defaultSigmaD013[ipt];
239  }
240  else {
241  AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
242  }
243  break;
244  default :
245  Warning("AliAnalysisTaskSEDmesonsFilterCJ", "%d not accepted!!", fCandidateType);
246  break;
247  }
248 
251 
252  DefineOutput(1, TList::Class()); // histos
253  DefineOutput(2, AliRDHFCuts::Class()); // my cuts
254  DefineOutput(3, TClonesArray::Class()); //array of candidates
255  DefineOutput(4, TClonesArray::Class()); //array of SB candidates
256  DefineOutput(5, TClonesArray::Class()); //array of candidates and event tracks
257  DefineOutput(6, TClonesArray::Class()); //array of SB candidates and event tracks
258  DefineOutput(7, TClonesArray::Class()); //array of MC D and event MC particles
259 }
260 
261 //_______________________________________________________________________________
263 {
264  //
265  // Destructor
266  //
267 
268  Info("~AliAnalysisTaskSEDmesonsFilterCJ","Calling Destructor");
269 
270  if (fCuts) {
271  delete fCuts;
272  fCuts = 0;
273  }
274 
275  if (fCandidateArray) {
276  delete fCandidateArray;
277  fCandidateArray = 0;
278  }
279 
280  if (fSideBandArray) {
281  delete fSideBandArray;
282  fSideBandArray = 0;
283  }
284 
285  if (fCombinedDmesons) {
286  delete fCombinedDmesons;
287  fCombinedDmesons = 0;
288  delete fMCCombinedDmesons;
289  fMCCombinedDmesons = 0;
290  }
291 
292  if (fCombinedDmesonsBkg) {
293  delete fCombinedDmesonsBkg;
295  }
296 }
297 
298 //_______________________________________________________________________________
300 {
301  //
302  // Initialization
303  //
304 
305  Info("AnalysisTaskSEDmesonsForJetCorrelations::Init()", "Entering method");
306 
307  switch (fCandidateType) {
308  case 0:
309  {
310  AliRDHFCutsD0toKpi* copyfCutsDzero = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
311  copyfCutsDzero->SetName("AnalysisCutsDzero");
312  PostData(2, copyfCutsDzero); // Post the data
313  } break;
314  case 1:
315  {
316  AliRDHFCutsDStartoKpipi* copyfCutsDstar = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
317  copyfCutsDstar->SetName("AnalysisCutsDStar");
318  PostData(2, copyfCutsDstar); // Post the cuts
319  } break;
320  default:
321  return;
322  }
323 
324  return;
325 }
326 
327 //_______________________________________________________________________________
329 {
330  //
331  // Create output objects
332  //
333 
334  Info("UserCreateOutputObjects","CreateOutputObjects of task %s", GetName());
335 
337 
338  DefineHistoForAnalysis(); // define histograms
339 
340  if (fUseReco) {
341  if (fCandidateType == kD0toKpi){
342  fCandidateArray = new TClonesArray("AliAODRecoDecayHF2Prong",10);
343  fSideBandArray = new TClonesArray("AliAODRecoDecayHF2Prong",10);
344  }
345  else if (fCandidateType == kDstartoKpipi) {
346  fCandidateArray = new TClonesArray("AliAODRecoCascadeHF",10);
347  fSideBandArray = new TClonesArray("AliAODRecoCascadeHF",10);
348  }
349  else {
350  AliWarning(Form("Candidate type %d not recognized!", fCandidateType));
351  return;
352  }
353  }
354  else {
355  fCandidateArray = new TClonesArray("AliAODMCParticle",10);
356  fSideBandArray = new TClonesArray("TObject",0); // not used
357  }
358 
359  if (fCombineDmesons) {
360  fCombinedDmesons = new TClonesArray("AliEmcalParticle",50);
361  fCombinedDmesonsBkg = new TClonesArray("AliEmcalParticle",50);
362  fMCCombinedDmesons = new TClonesArray("AliAODMCParticle",50);
363  }
364  else {
365  fCombinedDmesons = new TClonesArray("TObject",0); // not used
366  fCombinedDmesonsBkg = new TClonesArray("TObject",0); // not used
367  fMCCombinedDmesons = new TClonesArray("TObject",0); // not used
368  }
369 
370  fCandidateArray->SetOwner();
371  fCandidateArray->SetName(GetOutputSlot(3)->GetContainer()->GetName());
372 
373  //this is used for the DStar side bands and MC!
374  fSideBandArray->SetOwner();
375  fSideBandArray->SetName(GetOutputSlot(4)->GetContainer()->GetName());
376 
377  fCombinedDmesons->SetOwner();
378  fCombinedDmesons->SetName(GetOutputSlot(5)->GetContainer()->GetName());
379 
380  fCombinedDmesonsBkg->SetOwner();
381  fCombinedDmesonsBkg->SetName(GetOutputSlot(6)->GetContainer()->GetName());
382 
383  fMCCombinedDmesons->SetOwner();
384  fMCCombinedDmesons->SetName(GetOutputSlot(7)->GetContainer()->GetName());
385 
386  PostData(1, fOutput);
387  PostData(3, fCandidateArray);
388  PostData(4, fSideBandArray);
389  PostData(5, fCombinedDmesons);
390  PostData(6, fCombinedDmesonsBkg);
391  PostData(7, fMCCombinedDmesons);
392 
393  Info("UserCreateOutputObjects","Data posted for task %s", GetName());
394 }
395 
396 //_______________________________________________________________________________
398 {
399  //
400  // To be executed only once, for the first event
401  //
402 
403  AliDebug(2, "Entering ExecOnce()");
404 
405  if (fInhibitTask) return;
406 
407  // Load the event
408  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
409 
410  if (fAodEvent) {
411  fArrayDStartoD0pi = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(fBranchName.Data()));
412  }
413  else {
414  if (AODEvent() && IsStandardAOD()) {
415 
416  // In case there is an AOD handler writing a standard AOD, use the AOD
417  // event in memory rather than the input (ESD) event.
418  fAodEvent = dynamic_cast<AliAODEvent*>(AODEvent());
419 
420  // in this case the branches in the deltaAOD (AliAOD.VertexingHF.root)
421  // have to taken from the AOD event hold by the AliAODExtension
422  AliAODHandler *aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
423  if(aodHandler->GetExtensions()) {
424  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
425  AliAODEvent *aodFromExt = ext->GetAOD();
426  fArrayDStartoD0pi = (TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
427  }
428  else {
429  AliError(Form("This task need an AOD event! Task '%s' will be disabled!", GetName()));
430  fInhibitTask = kTRUE;
431  return;
432  }
433  }
434  }
435 
436  if (fArrayDStartoD0pi) {
437  TString objname(fArrayDStartoD0pi->GetClass()->GetName());
438  TClass cls(objname);
439  if (!cls.InheritsFrom("AliAODRecoDecayHF2Prong")) {
440  AliError(Form("%s: Objects of type %s in %s are not inherited from AliAODRecoDecayHF2Prong! Task will be disabled!",
441  GetName(), cls.GetName(), fArrayDStartoD0pi->GetName()));
442  fInhibitTask = kTRUE;
443  fArrayDStartoD0pi = 0;
444  return;
445  }
446  }
447  else {
448  AliError(Form("Could not find array %s, skipping the event. Task '%s' will be disabled!", fBranchName.Data(), GetName()));
449  fInhibitTask = kTRUE;
450  return;
451  }
452 
453  if (fUseMCInfo) {
454  fMCarray = dynamic_cast<TClonesArray*>(fAodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
455  if (!fMCarray) {
456  AliError(Form("MC particles not found! Task '%s' will be disabled!", GetName()));
457  fInhibitTask = kTRUE;
458  return;
459  }
460 
461  fMCHeader = (AliAODMCHeader*)fAodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
462  if(!fMCHeader) {
463  AliError(Form("MC header not found! Task '%s' will be disabled!", GetName()));
464  return;
465  }
466  }
467 
468  if (fCombineDmesons) {
472  }
473 
474  fRan = new TRandom3(0);
475 
477 }
478 
479 //_______________________________________________________________________________
481 {
482  //
483  // Analysis execution
484  //
485 
486  if (fInhibitTask) return kFALSE;
487 
488  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
489  if (matchingAODdeltaAODlevel<=0) {
490  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
491  return kFALSE;
492  }
493 
494  AliDebug(2, "Entering Run()");
495 
496  //clear the TClonesArray from the previous event
497  fCandidateArray->Clear();
498  fSideBandArray->Clear();
499  fCombinedDmesons->Clear();
500  fCombinedDmesonsBkg->Clear();
501  AliDebug(2, "TClonesArray cleared");
502 
503  fHistStat->Fill(0);
504 
505  // fix for temporary bug in ESDfilter
506  // the AODs with null vertex pointer didn't pass the PhysSel
507  if (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001) return kFALSE;
508 
509  //Event selection
510  Bool_t iseventselected = fCuts->IsEventSelected(fAodEvent);
511  if (!iseventselected) return kFALSE;
512  fHistStat->Fill(1);
513 
514  AliDebug(2, "Event selected");
515 
516  const Int_t nD = fArrayDStartoD0pi->GetEntriesFast();
517  AliDebug(2, Form("Found %d vertices", nD));
518  if (!fUseMCInfo) fHistStat->Fill(2, nD);
519 
520 
521  Int_t pdgMeson = 413;
522  if (fCandidateType == kD0toKpi) pdgMeson = 421;
523 
524  fNCand = 0;
525  fNSBCand = 0;
526 
528 
529 
530 // for MC response matrix of efficiency studies, fMultCand option only
531 if(fUseMCInfo && fBuildRMEff){
532 
533  const Int_t nDMC = fMCarray->GetEntriesFast();
534 
535  for (Int_t iMCcharm = 0; iMCcharm < nDMC; iMCcharm++){ //loop over MC D
536 
537  AliAODMCParticle* charmPart = static_cast<AliAODMCParticle*>(fMCarray->At(iMCcharm));
538  if(TMath::Abs(charmPart->GetPdgCode()) != pdgMeson) continue;
539  if(!charmPart) continue;
540  fHistStat->Fill(2);
541 
542  Int_t origin = CheckOrigin(charmPart, fMCarray);
543  if (origin < 0) continue;
544  if (fRejectQuarkNotFound && origin == kQuarkNotFound) {
545  fHistStat->Fill(6);
546  continue;
547  }
548  if (fRejectDfromB && origin == kFromBottom) {
549  fHistStat->Fill(7);
550  continue;
551  }
552  if (fKeepOnlyDfromB && origin != kFromBottom) {
553  fHistStat->Fill(8);
554  continue;
555  }
556 
557  Int_t decay = CheckDecayChannel(charmPart, fMCarray);
558  if( TMath::Abs(charmPart->GetPdgCode()) == 413 && decay != kDecayDStartoKpipi) continue;
559  if( TMath::Abs(charmPart->GetPdgCode()) == 421 && decay != kDecayD0toKpi) continue;
560 
561  fHistStat->Fill(3);
562 
563  if (fNCand==fAnalyseCand){
564 
565  new ((*fMCCombinedDmesons)[0]) AliAODMCParticle(*charmPart);
566 
567  AliAODRecoDecayHF2Prong* charmCand = 0;
568  AliAODRecoCascadeHF* dstar = 0;
569 
570  // loop over reco D candidates to find a match to MC
571  Int_t isRecoD = kFALSE;
572  for (Int_t icharm = 0; icharm < nD; icharm++) {
573 
574  charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fArrayDStartoD0pi->At(icharm)); // D candidates
575  if (!charmCand) continue;
576  if(!(vHF->FillRecoCand(fAodEvent,charmCand))) continue;
577 
578  Int_t nprongs = charmCand->GetNProngs();
579  AliDebug(2, Form("Candidate is %d, and nprongs = %d", fCandidateType, nprongs));
580 
581  if (fCandidateType == kDstartoKpipi) {
582  dstar = dynamic_cast<AliAODRecoCascadeHF*>(charmCand);
583  if (!dstar) {
584  Error("AliAnalysisTaskSEDmesonsFilterCJ::UserExec","Candidate type is D* but object type is wrong (should be AliAODRecoCascadeHF)");
585  continue;
586  }
587  }
588 
589  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
590  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
591 
592  Int_t mcLabel = 0;
593  if (fCandidateType == kDstartoKpipi) mcLabel = dstar->MatchToMC(413, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCarray);
594  else mcLabel = charmCand->MatchToMC(421, fMCarray, fNProngs, fPDGdaughters);
595 
596  if(mcLabel == iMCcharm) { isRecoD = kTRUE; break; }
597  }
598 
599  if (!isRecoD) break;
600  if (!charmCand) break;
601  if (fCandidateType == kDstartoKpipi && !dstar) break;
602 
603  fHistStat->Fill(4);
604 
605  // region of interest + cuts
606  if (!fCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(pdgMeson))) break;
607 
608  //candidate selected by cuts and PID
609  Int_t isSelected = 0;
610  isSelected = fCuts->IsSelected(charmCand, AliRDHFCuts::kAll, fAodEvent); //selected
611  if (!isSelected) break;
612 
613  fHistStat->Fill(5);
614 
615  if (fCandidateType == kDstartoKpipi) {
616  AliAODRecoDecayHF2Prong* D0fromDstar = dstar->Get2Prong();
617  fHistInvMassS->Fill(dstar->DeltaInvMass());
618  fHistImpParS->Fill(dstar->Getd0Prong(0), dstar->PtProng(0)); //bachelor
619  fHistImpParS->Fill(D0fromDstar->Getd0Prong(0), D0fromDstar->PtProng(0));
620  fHistImpParS->Fill(D0fromDstar->Getd0Prong(1), D0fromDstar->PtProng(1));
621  }
622  else {
623  fHistImpParS->Fill(charmCand->Getd0Prong(0), charmCand->PtProng(0));
624  fHistImpParS->Fill(charmCand->Getd0Prong(1), charmCand->PtProng(1));
625  fHistInvMassS->Fill(charmCand->InvMassD0());
626  fHistInvMassS->Fill(charmCand->InvMassD0bar());
627  }
628 
629  if (fCandidateType == kDstartoKpipi) {
630  Int_t isDstar = charmPart == NULL ? 0 : 1;
631  //fill histograms of kinematics, using MC truth
632  FillDStarMCTruthKinHistos(dstar, isSelected, isDstar);
633 
634  new ((*fCandidateArray)[0]) AliAODRecoCascadeHF(*dstar);
635  new ((*fCombinedDmesons)[0]) AliEmcalParticle(dstar);
636 
637  }
638  else {
639 
640  Int_t isD0 = 0;
641  Int_t pdgCode = charmPart->GetPdgCode();
642  if (pdgCode == 421) { isD0 = 1; }
643  else if (pdgCode == -421) { isD0 = -1; }
644  //fill histograms of kinematics, using MC truth
645  FillD0MCTruthKinHistos(charmCand, isSelected, isD0);
646 
647  new ((*fCandidateArray)[0]) AliAODRecoDecayHF2Prong(*charmCand);
648  new ((*fCombinedDmesons)[0]) AliEmcalParticle(charmCand);
649  }
650 
651  break;
652  }
653  else { fNCand++; }
654 
655 }
656 }
657 // for data, or MC without RM or efficiency studies
658 else {
659  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
660  Int_t isSelected = 0;
661 
662  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fArrayDStartoD0pi->At(icharm)); // D candidates
663  if (!charmCand) continue;
664  if(!(vHF->FillRecoCand(fAodEvent,charmCand))) continue;
665 
666  Int_t nprongs = charmCand->GetNProngs();
667  AliDebug(2, Form("Candidate is %d, and nprongs = %d", fCandidateType, nprongs));
668 
669  AliAODRecoCascadeHF* dstar = 0;
670 
671  if (fCandidateType == kDstartoKpipi) {
672  dstar = dynamic_cast<AliAODRecoCascadeHF*>(charmCand);
673  if (!dstar) {
674  Error("AliAnalysisTaskSEDmesonsFilterCJ::UserExec","Candidate type is D* but object type is wrong (should be AliAODRecoCascadeHF)");
675  continue;
676  }
677  }
678 
679  // region of interest + cuts
680  if (!fCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(pdgMeson))) continue;
681 
683  FillDstarSideBands(dstar);
684  }
685 
686 
687  //candidate selected by cuts and PID
688  isSelected = fCuts->IsSelected(charmCand, AliRDHFCuts::kAll, fAodEvent); //selected
689  if (!isSelected) continue;
690 
691  if (fCandidateType == kDstartoKpipi) {
692  ProcessDstar(dstar, isSelected);
693  }
694  else {
695  ProcessD0(charmCand, isSelected);
696  }
697  } // end of D cand loop
698 }
699 
700 
701  delete vHF;
702 
703 
704  AliDebug(2, "Loop done");
705 
706  if (fCombineDmesons) {
707  if (fCombinedDmesons->GetEntriesFast() > 0) {
709  }
710 
711  if (fMCCombinedDmesons->GetEntriesFast() > 0) {
713  }
714 
715  if (fCombinedDmesonsBkg->GetEntriesFast() > 0) {
717  }
718  }
719 
720  fHistNCandEv->Fill(fCandidateArray->GetEntriesFast());
722  Int_t nsbcand = fSideBandArray->GetEntriesFast();
723  fHistStat->Fill(4, nsbcand);
724  fHistNSBCandEv->Fill(nsbcand);
725  }
726 
727  PostData(1, fOutput);
728  PostData(3, fCandidateArray);
729  PostData(4, fSideBandArray);
730  PostData(5, fCombinedDmesons);
731  PostData(6, fCombinedDmesonsBkg);
732  PostData(7, fMCCombinedDmesons);
733 
734  AliDebug(2, "Exiting method");
735 
736  return kTRUE;
737 }
738 
739 //_______________________________________________________________________________
741 {
742  // Process the D0 candidate.
743 
744  // For MC analysis look for the AOD MC particle associated with the charm candidate
745  AliAODMCParticle* charmPart = 0;
746  if (fUseMCInfo) {
747  Int_t mcLabel = charmCand->MatchToMC(421, fMCarray, fNProngs, fPDGdaughters);
748 
749  if (mcLabel >= 0) {
750  charmPart = static_cast<AliAODMCParticle*>(fMCarray->At(mcLabel));
751  }
752  }
753 
754  if (charmPart) {
755 
756  Int_t origin = CheckOrigin(charmPart, fMCarray);
757  if (origin < 0) return;
758 
759  if (fRejectQuarkNotFound && origin == kQuarkNotFound) {
760  fHistStat->Fill(6);
761  return;
762  }
763  if (fRejectDfromB && origin == kFromBottom) {
764  fHistStat->Fill(7);
765  return;
766  }
767  if (fKeepOnlyDfromB && origin != kFromBottom) {
768  fHistStat->Fill(8);
769  return;
770  }
771 
772  fHistStat->Fill(2);
773  }
774  else {
775  fHistStat->Fill(5);
776  }
777 
778  // For MC background fill fSideBandArray
779  if (fUseMCInfo && !charmPart) {
780  if (fUseReco) {
781  if(!fMultCand) new ((*fSideBandArray)[fNSBCand]) AliAODRecoDecayHF2Prong(*charmCand);
782  else if(fNSBCand==fAnalyseCand) new ((*fSideBandArray)[0]) AliAODRecoDecayHF2Prong(*charmCand);
783 
784  if (fCombineDmesons) {
785  if(!fMultCand) new ((*fCombinedDmesonsBkg)[fNSBCand]) AliEmcalParticle(charmCand);
786  else if(fNSBCand==fAnalyseCand) new ((*fCombinedDmesonsBkg)[0]) AliEmcalParticle(charmCand);
787  }
788 
789  fHistImpParB->Fill(charmCand->Getd0Prong(0), charmCand->PtProng(0));
790  fHistImpParB->Fill(charmCand->Getd0Prong(1), charmCand->PtProng(1));
791  fHistInvMassB->Fill(charmCand->InvMassD0());
792  fHistInvMassB->Fill(charmCand->InvMassD0bar());
793 
794  fNSBCand++;
795  }
796  }
797  // For data or MC signal fill fCandidateArray
798  else {
799  // For data or MC with the requirement fUseReco fill with candidates
800  if (fUseReco) {
801  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODRecoDecayHF2Prong(*charmCand);
802  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODRecoDecayHF2Prong(*charmCand);
803 
804  if (fCombineDmesons) {
805  if(!fMultCand) new ((*fCombinedDmesons)[fNCand]) AliEmcalParticle(charmCand);
806  else if(fNCand==fAnalyseCand)
807  {
808  new ((*fCombinedDmesons)[0]) AliEmcalParticle(charmCand);
809  if(charmPart) new ((*fMCCombinedDmesons)[0]) AliAODMCParticle(*charmPart);
810  }
811  }
812 
813  fHistImpParS->Fill(charmCand->Getd0Prong(0), charmCand->PtProng(0));
814  fHistImpParS->Fill(charmCand->Getd0Prong(1), charmCand->PtProng(1));
815  fHistInvMassS->Fill(charmCand->InvMassD0());
816  fHistInvMassS->Fill(charmCand->InvMassD0bar());
817  }
818  // For MC with requirement particle level fill with AliAODMCParticle
819  else {
820  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODMCParticle(*charmPart);
821  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODMCParticle(*charmPart);
822  }
823  fHistStat->Fill(3);
824  fNCand++;
825  }
826 
827  // Now filling some more histograms
828 
829  // mass vs pt
830  if (isSelected == 1 || isSelected == 3) fHistInvMassPtD->Fill(charmCand->InvMassD0(), charmCand->Pt());
831  if (isSelected >= 2) fHistInvMassPtD->Fill(charmCand->InvMassD0bar(), charmCand->Pt());
832 
833  if (fUseMCInfo) { //fill histograms of kinematics, using MC truth
834  Int_t isD0 = 0;
835  if (charmPart) {
836  Int_t pdgCode = charmPart->GetPdgCode();
837 
838  if (pdgCode == 421) {
839  isD0 = 1;
840  }
841  else if (pdgCode == -421) {
842  isD0 = -1;
843  }
844  else {
845  AliDebug(2, "Not a D0/D0bar meson!");
846  return;
847  }
848  }
849 
850  FillD0MCTruthKinHistos(charmCand, isSelected, isD0);
851  }
852 }
853 
854 //_______________________________________________________________________________
856 {
857  // Process the D* candidate.
858 
859  AliDebug(2, "Entering method");
860 
861  // For MC analysis look for the AOD MC particle associated with the charm candidate
862  AliAODMCParticle* charmPart = 0;
863  if (fUseMCInfo) {
864  //D* and D0 prongs needed to MatchToMC method
865  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
866  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
867 
868  Int_t mcLabel = dstar->MatchToMC(413, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCarray);
869 
870  if (mcLabel >= 0) {
871  charmPart = static_cast<AliAODMCParticle*>(fMCarray->At(mcLabel));
872  }
873  }
874 
875  if (charmPart) {
876 
877  Int_t origin = CheckOrigin(charmPart, fMCarray);
878  if (origin < 0) return;
879 
880  if (fRejectQuarkNotFound && origin == kQuarkNotFound) {
881  fHistStat->Fill(6);
882  return;
883  }
884  if (fRejectDfromB && origin == kFromBottom) {
885  fHistStat->Fill(7);
886  return;
887  }
888  if (fKeepOnlyDfromB && origin != kFromBottom) {
889  fHistStat->Fill(8);
890  return;
891  }
892 
893  fHistStat->Fill(2);
894  }
895  else {
896  fHistStat->Fill(5);
897  }
898 
899  AliAODRecoDecayHF2Prong* D0fromDstar = dstar->Get2Prong();
900 
901  // For MC background fill fSideBandArray
902  if (fUseMCInfo && !charmPart) {
903  if (fUseReco) {
904  if(!fMultCand) new ((*fSideBandArray)[fNSBCand]) AliAODRecoCascadeHF(*dstar);
905  else if(fNSBCand==fAnalyseCand) new ((*fSideBandArray)[0]) AliAODRecoCascadeHF(*dstar);
906 
907  if (fCombineDmesons) {
908  if(!fMultCand) new ((*fCombinedDmesonsBkg)[fNSBCand]) AliEmcalParticle(dstar);
909  else if(fNSBCand==fAnalyseCand) new ((*fCombinedDmesonsBkg)[0]) AliEmcalParticle(dstar);
910  }
911 
912  fHistInvMassB->Fill(dstar->DeltaInvMass());
913  fHistImpParB->Fill(dstar->Getd0Prong(0), dstar->PtProng(0)); //bachelor
914  fHistImpParB->Fill(D0fromDstar->Getd0Prong(0), D0fromDstar->PtProng(0));
915  fHistImpParB->Fill(D0fromDstar->Getd0Prong(1), D0fromDstar->PtProng(1));
916 
917  fNSBCand++;
918  }
919  }
920  // For data and MC signal fill fCandidateArray
921  else {
922  // For data or MC signal with the requirement fUseReco fill with candidates
923  if (fUseReco) {
924  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODRecoCascadeHF(*dstar);
925  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODRecoCascadeHF(*dstar);
926 
927  if (fCombineDmesons) {
928  if(!fMultCand) new ((*fCombinedDmesons)[fNCand]) AliEmcalParticle(dstar);
929  else if(fNCand==fAnalyseCand)
930  {
931  new ((*fCombinedDmesons)[0]) AliEmcalParticle(dstar);
932  if(charmPart) new ((*fMCCombinedDmesons)[0]) AliAODMCParticle(*charmPart);
933  }
934  }
935  }
936  // For MC signal with requirement particle level fill with AliAODMCParticle
937  else {
938  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODMCParticle(*charmPart);
939  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODMCParticle(*charmPart);
940  }
941  fNCand++;
942 
943  fHistInvMassS->Fill(dstar->DeltaInvMass());
944  fHistImpParS->Fill(dstar->Getd0Prong(0), dstar->PtProng(0)); //bachelor
945  fHistImpParS->Fill(D0fromDstar->Getd0Prong(0), D0fromDstar->PtProng(0));
946  fHistImpParS->Fill(D0fromDstar->Getd0Prong(1), D0fromDstar->PtProng(1));
947  fHistStat->Fill(3);
948  }
949 
950  // Now filling some more histograms.
951 
952  // select D* in the D0 window.
953  // In the cut object window is loose to allow for side bands
954 
955  // retrieve the corresponding pt bin for the candidate
956  // and set the expected D0 width (x3)
957 
958  Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
959  Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
960  Double_t invMassPDG = mPDGDstar - mPDGD0;
961 
962  Double_t ptD = dstar->Pt();
963  Int_t ptDbin = fCuts->PtBin(ptD);
964  if (ptDbin < 0 || ptDbin >= fCuts->GetNPtBins()) {
965  AliError(Form("Pt %.3f out of bounds", ptD));
966  return;
967  }
968 
969  AliDebug(1, Form("Pt bin %d and sigma D0 %.4f", ptDbin, fSigmaD0[ptDbin]));
970  //consider the D* candidates only if the mass of the D0 is in 3 sigma wrt the PDG value
971  if ((dstar->InvMassD0()>=(mPDGD0-3.*fSigmaD0[ptDbin])) && (dstar->InvMassD0()<=(mPDGD0+3.*fSigmaD0[ptDbin]))) {
972  AliDebug(2, "D0 mass within 3 sigma");
973 
974  // D* delta mass
975  AliDebug(2, "Filling invariant mass vs pt histogram");
976  fHistInvMassPtD->Fill(dstar->DeltaInvMass(), ptD); // 2 D slice for pt bins
977 
978  // Soft pion pt for good candidates
979  Double_t invmassDelta = dstar->DeltaInvMass();
980  if (TMath::Abs(invmassDelta - invMassPDG) < 0.0021) {
981  AliDebug(2, "Filling pion pt histogram");
982  //softpion from D* decay
983  AliAODTrack *softPionTrack = static_cast<AliAODTrack*>(dstar->GetBachelor());
984  if (softPionTrack) fHistPtPion->Fill(softPionTrack->Pt());
985  }
986  }
987 
988  if (fUseMCInfo) { //fill histograms of kinematics, using MC truth
989  Int_t isDstar = charmPart == NULL ? 0 : 1;
990  FillDStarMCTruthKinHistos(dstar, isDstar, isSelected);
991  }
992 
993  AliDebug(2, "Exiting method");
994 }
995 
996 //_______________________________________________________________________________
998 {
999  // Fill some histogram with kinematic information of the D0 candidate.
1000 
1001  Double_t ptD = charmCand->Pt();
1002  Double_t aD = charmCand->Phi();
1003  Double_t adaugh[2] = {charmCand->PhiProng(0), charmCand->PhiProng(1)};
1004  AliAODTrack* p0 = static_cast<AliAODTrack*>(charmCand->GetDaughter(0));
1005  AliAODTrack* p1 = static_cast<AliAODTrack*>(charmCand->GetDaughter(1));
1006  Float_t dR0 = DeltaR(charmCand, p0);
1007  Float_t dR1 = DeltaR(charmCand, p1);
1008 
1009  if (isD0 == 0) { //background
1010  if (isSelected == 1 || isSelected == 3) { // selected as D0
1011  fHistAlphaDKB->Fill(aD-adaugh[0],ptD);
1012  fHistAlphaDpiB->Fill(aD-adaugh[1],ptD);
1013 
1014  fHistDeltaRDKB->Fill(dR0,ptD);
1015  fHistDeltaRDpiB->Fill(dR1,ptD);
1016  }
1017  if (isSelected >= 2) { //selected as D0bar
1018  fHistAlphaDpiB->Fill(aD-adaugh[0],ptD);
1019  fHistAlphaDKB->Fill(aD-adaugh[1],ptD);
1020 
1021  fHistDeltaRDpiB->Fill(dR0,ptD);
1022  fHistDeltaRDKB->Fill(dR1,ptD);
1023  }
1024  }
1025  else if (isD0 == 1) { //D0
1026  fHistAlphaDKS->Fill(aD-adaugh[0],ptD);
1027  fHistAlphaDpiS->Fill(aD-adaugh[1],ptD);
1028 
1029  fHistDeltaRDKS->Fill(dR0,ptD);
1030  fHistDeltaRDpiS->Fill(dR1,ptD);
1031 
1032  if (isSelected == 3) { // selected as both D0bar/D0
1033  fHistAlphaDpiR->Fill(aD-adaugh[0],ptD);
1034  fHistAlphaDKR->Fill(aD-adaugh[1],ptD);
1035 
1036  fHistDeltaRDpiR->Fill(dR0,ptD);
1037  fHistDeltaRDKR->Fill(dR1,ptD);
1038  }
1039  }
1040  else if (isD0 == -1) { //D0bar
1041  fHistAlphaDKS->Fill(aD-adaugh[1],ptD);
1042  fHistAlphaDpiS->Fill(aD-adaugh[0],ptD);
1043 
1044  fHistDeltaRDKS->Fill(dR1,ptD);
1045  fHistDeltaRDpiS->Fill(dR0,ptD);
1046 
1047  if (isSelected == 3) { // selected as both D0bar/D0
1048  fHistAlphaDpiR->Fill(aD-adaugh[1],ptD);
1049  fHistAlphaDKR->Fill(aD-adaugh[0],ptD);
1050 
1051  fHistDeltaRDpiR->Fill(dR1, ptD);
1052  fHistDeltaRDKR->Fill(dR0, ptD);
1053  }
1054  }
1055 }
1056 
1057 //_______________________________________________________________________________
1059 {
1060  // Fill some histogram with kinematic information of the D0 candidate.
1061 
1062  AliAODTrack *softPionTrack = static_cast<AliAODTrack*>(dstar->GetBachelor());
1063  Double_t ptD = dstar->Pt();
1064  Double_t aD = dstar->Phi();
1065  Double_t apis= softPionTrack->Phi();
1066 
1067  AliAODRecoDecayHF2Prong* D0fromDstar = dstar->Get2Prong();
1068  Double_t aD0 = D0fromDstar->Phi();
1069  Int_t isD0 = D0fromDstar->Charge()>0 ? kTRUE : kFALSE;
1070  Double_t aK = isD0 ? D0fromDstar->PhiProng(0) : D0fromDstar->PhiProng(1);
1071  Double_t api = isD0 ? D0fromDstar->PhiProng(1) : D0fromDstar->PhiProng(0);
1072  Double_t dRDD0 = DeltaR(dstar,D0fromDstar);
1073  Double_t dRDpis = DeltaR(dstar,softPionTrack);
1074  Double_t dRDpi = DeltaR(dstar, isD0 ? static_cast<AliVParticle*>(D0fromDstar->GetDaughter(1)) : static_cast<AliVParticle*>(D0fromDstar->GetDaughter(0)));
1075  Double_t dRDK = DeltaR(dstar, isD0 ? static_cast<AliVParticle*>(D0fromDstar->GetDaughter(0)) : static_cast<AliVParticle*>(D0fromDstar->GetDaughter(1)));
1076 
1077  if (isDstar) {
1078  fHistAlphaDDS ->Fill(aD-aD0, ptD);
1079  fHistAlphaDpisS->Fill(aD-apis, ptD);
1080  fHistAlphaDpiS ->Fill(aD-api, ptD);
1081  fHistAlphaDKS ->Fill(aD-aK, ptD);
1082 
1083  fHistDeltaRDDS ->Fill(dRDD0, ptD);
1084  fHistDeltaRDpisS->Fill(dRDpis, ptD);
1085  fHistDeltaRDpiS ->Fill(dRDpi, ptD);
1086  fHistDeltaRDKS ->Fill(dRDK, ptD);
1087  }
1088  else {
1089  fHistAlphaDDB ->Fill(aD-aD0, ptD);
1090  fHistAlphaDpisB->Fill(aD-apis, ptD);
1091  fHistAlphaDpiB ->Fill(aD-api, ptD);
1092  fHistAlphaDKB ->Fill(aD-aK, ptD);
1093 
1094  fHistDeltaRDDB ->Fill(dRDD0, ptD);
1095  fHistDeltaRDpisB->Fill(dRDpis, ptD);
1096  fHistDeltaRDpiB ->Fill(dRDpi, ptD);
1097  fHistDeltaRDKB ->Fill(dRDK, ptD);
1098  }
1099 }
1100 
1101 //_______________________________________________________________________________
1103 {
1104  // Fills the array of Dstar side band candidates.
1105 
1106  //select by track cuts the side band candidates (don't want mass cut)
1107  Int_t isSelected = fCuts->IsSelected(dstar, AliRDHFCuts::kTracks, fAodEvent);
1108  if (!isSelected) return;
1109 
1110  //add a reasonable cut on the invariant mass (e.g. (+-2\sigma, +-10 \sigma), with \sigma = fSigmaD0[bin])
1111  Int_t bin = fCuts->PtBin(dstar->Pt());
1112  if (bin < 0 || bin >= fCuts->GetNPtBins()) {
1113  AliError(Form("Pt %.3f out of bounds", dstar->Pt()));
1114  return;
1115  }
1116 
1117  const Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1118 
1119  //if data and Dstar from D0 side band
1120  if (((dstar->InvMassD0() < (mPDGD0-3.*fSigmaD0[bin])) && (dstar->InvMassD0() > (mPDGD0-10.*fSigmaD0[bin]))) /*left side band*/ ||
1121  ((dstar->InvMassD0() > (mPDGD0+3.*fSigmaD0[bin])) && (dstar->InvMassD0() < (mPDGD0+10.*fSigmaD0[bin]))) /*right side band*/) {
1122 
1123  if(!fMultCand) new ((*fSideBandArray)[fNSBCand]) AliAODRecoCascadeHF(*dstar);
1124  else if(fNSBCand==fAnalyseCand) new ((*fSideBandArray)[0]) AliAODRecoCascadeHF(*dstar);
1125  fNSBCand++;
1126  }
1127 }
1128 
1129 //_______________________________________________________________________________
1131 {
1132  // Set the mass limits using the PDG code and the given range.
1133 
1134  Float_t mass = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
1135 
1136  // compute the Delta mass for the D*
1137  if (fCandidateType==kDstartoKpipi) mass -= TDatabasePDG::Instance()->GetParticle(421)->Mass();
1138 
1139  fMinMass = mass - range;
1140  fMaxMass = mass + range;
1141 
1142  AliInfo(Form("Setting mass limits to %f, %f", fMinMass, fMaxMass));
1143  if ((fMinMass<0.) || (fMaxMass<=0.) || (fMaxMass<fMinMass)) {
1144  AliError(Form("Wrong mass limits! Task '%s' will be disabled!", GetName()));
1145  fInhibitTask = kTRUE;
1146  }
1147 }
1148 
1149 //_______________________________________________________________________________
1151 {
1152  // Set the mass limits.
1153 
1154  if (uplimit>lowlimit) {
1155  fMinMass = lowlimit;
1156  fMaxMass = uplimit;
1157  }
1158  else {
1159  printf("Error! Lower limit larger than upper limit!\n");
1160  fMinMass = uplimit - uplimit*0.2;
1161  fMaxMass = uplimit;
1162  }
1163 }
1164 
1165 //_______________________________________________________________________________
1167 {
1168  // Set the D0 width for the D* analysis.
1169 
1170  if (nptbins > 30) {
1171  AliWarning("Maximum number of bins allowed is 30!");
1172  return kFALSE;
1173  }
1174 
1175  if (!width) return kFALSE;
1176  for (Int_t ipt=0; ipt < nptbins; ipt++) fSigmaD0[ipt] = width[ipt];
1177 
1178  return kTRUE;
1179 }
1180 
1181 //_______________________________________________________________________________
1183 {
1184  // Allocate the histograms for the analysis.
1185 
1186  // Statistics
1187  fHistStat = new TH1I("fHistStat", "Statistics", 9, -0.5, 8.5);
1188  fHistStat->GetXaxis()->SetBinLabel(1, "N ev anal");
1189  fHistStat->GetXaxis()->SetBinLabel(2, "N ev sel");
1190  if(fUseMCInfo && fBuildRMEff){
1191  fHistStat->GetXaxis()->SetBinLabel(3, "N Gen D");
1192  fHistStat->GetXaxis()->SetBinLabel(4, "N Gen Sel D");
1193  fHistStat->GetXaxis()->SetBinLabel(5, "N D");
1194  fHistStat->GetXaxis()->SetBinLabel(6, "N cand sel cuts");
1195  fHistStat->GetXaxis()->SetBinLabel(7, "N rej no quark");
1196  fHistStat->GetXaxis()->SetBinLabel(8, "N rej from B");
1197  fHistStat->GetXaxis()->SetBinLabel(9, "N rej from D");
1198  }
1199  else {
1200  if (fUseMCInfo) {
1201  if (fUseReco) {
1202  fHistStat->GetXaxis()->SetBinLabel(3, "N D");
1203  }
1204  else {
1205  fHistStat->GetXaxis()->SetBinLabel(3, "N Gen D");
1206  }
1207  }
1208  else {
1209  fHistStat->GetXaxis()->SetBinLabel(3, "N Cand");
1210  }
1211  fHistStat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
1212  if (fCandidateType == kDstartoKpipi) {
1213  fHistStat->GetXaxis()->SetBinLabel(5, "N side band cand");
1214  }
1215  if (fUseMCInfo) {
1216  fHistStat->GetXaxis()->SetBinLabel(6, "N Background");
1217  fHistStat->GetXaxis()->SetBinLabel(7, "N rej no quark");
1218  fHistStat->GetXaxis()->SetBinLabel(8, "N rej from B");
1219  fHistStat->GetXaxis()->SetBinLabel(9, "N rej from D");
1220  }
1221  }
1222 
1223  fHistStat->SetNdivisions(1);
1224  fOutput->Add(fHistStat);
1225 
1226  fHistNCandEv = new TH1F("fHistNCandEv", "Number of candidates per event (after cuts);# cand/ev", 100, 0., 100.);
1227  fOutput->Add(fHistNCandEv);
1228 
1229  if ((fCandidateType == kDstartoKpipi) || fUseMCInfo) {
1230  fHistNSBCandEv = new TH1F("hnSBCandEv", "Number of side bands candidates per event (after cuts);# cand/ev", 100, 0.,100.);
1231  fOutput->Add(fHistNSBCandEv);
1232  }
1233 
1234  if (fCandidateType == kDstartoKpipi) {
1235  fHistPtPion = new TH1F("fHistPtPion", "Primary pions candidates pt", 500, 0., 10.);
1236  fHistPtPion->SetStats(kTRUE);
1237  fHistPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1238  fHistPtPion->GetYaxis()->SetTitle("entries");
1239  fOutput->Add(fHistPtPion);
1240  }
1241 
1242  // Invariant mass related histograms
1243  const Int_t nbinsmass = 200;
1244  const Int_t ptbinsD = 100;
1245  Float_t ptmin =0.;
1246  Float_t ptmax =50.;
1247  fHistInvMassPtD = new TH2F("fHistInvMassPtD", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, ptbinsD, ptmin, ptmax);
1248  fHistInvMassPtD->SetStats(kTRUE);
1249  fHistInvMassPtD->GetXaxis()->SetTitle("mass (GeV/c)");
1250  fHistInvMassPtD->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1251  fOutput->Add(fHistInvMassPtD);
1252 
1253  fHistImpParS = new TH2F("fHistImpParS", "Impact parameter of daughter tracks; Getd0Prong();#it{p}^{daugh}_{T} (GeV/c)",200, -0.1,0.1,ptbinsD, ptmin, ptmax); //same range of pt of the D, but pt daughter used
1254  fOutput->Add(fHistImpParS);
1255 
1256  fHistInvMassS = new TH1F("fHistInvMassS", "D invariant mass distribution (filled with fCandidateArray)", nbinsmass, fMinMass, fMaxMass);
1257  fHistInvMassS->SetStats(kTRUE);
1258  fHistInvMassS->GetXaxis()->SetTitle("mass (GeV/c)");
1259  fHistInvMassS->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1260  fOutput->Add(fHistInvMassS);
1261 
1262  const Int_t nbinsalpha = 200;
1263  Float_t minalpha = -TMath::Pi();
1264  Float_t maxalpha = TMath::Pi();
1265  const Int_t nbinsdeltaR = 200;
1266  Float_t mindeltaR = 0.;
1267  Float_t maxdeltaR = 10.;
1268 
1269 
1270 
1271  if (fUseMCInfo) {
1272  fHistImpParB = new TH2F("fHistImpParB", "Impact parameter of daughter tracks (Background); Getd0Prong();#it{p}^{daugh}_{T} (GeV/c)",200, -0.1,0.1,ptbinsD, ptmin, ptmax); //same range of pt of the D, but pt daughter used
1273  fOutput->Add(fHistImpParB);
1274 
1275  fHistInvMassB = new TH1F("fHistInvMassB", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass);
1276  fHistInvMassB->SetStats(kTRUE);
1277  fHistInvMassB->GetXaxis()->SetTitle("mass (GeV/c)");
1278  fHistInvMassB->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1279  fOutput->Add(fHistInvMassB);
1280 
1281  if (fCandidateType == kDstartoKpipi) {
1282 
1283  fHistAlphaDDS = new TH2F("fHistAlphaDDS","Angle D^{*}-D^{0} (Signal);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1284  fHistAlphaDpisS = new TH2F("fHistAlphaDpisS","Angle D^{*}-#pi_{soft} (Signal);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1285  fHistAlphaDpiS = new TH2F("fHistAlphaDpiS","Angle D^{*}-#pi (Signal);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1286  fHistAlphaDKS = new TH2F("fHistAlphaDKS","Angle D^{*}-K (Signal);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1287 
1288  fHistAlphaDDB = new TH2F("fHistAlphaDDB","Angle D^{*}-D^{0} (Background);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1289  fHistAlphaDpisB = new TH2F("fHistAlphaDpisB","Angle D^{*}-#pi_{soft} (Background);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1290  fHistAlphaDpiB = new TH2F("fHistAlphaDpiB","Angle D^{*}-#pi (Background);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1291  fHistAlphaDKB = new TH2F("fHistAlphaDKB","Angle D^{*}-K (Background);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1292 
1293  fHistDeltaRDDS = new TH2F("fHistDeltaRDDS","Angle D^{*}-D^{0} (Signal);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1294  fHistDeltaRDpisS = new TH2F("fHistDeltaRDpisS","Angle D^{*}-#pi_{soft} (Signal);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1295  fHistDeltaRDpiS = new TH2F("fHistDeltaRDpiS","Angle D^{*}-#pi (Signal);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1296  fHistDeltaRDKS = new TH2F("fHistDeltaRDKS","Angle D^{*}-K (Signal);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1297 
1298  fHistDeltaRDDB = new TH2F("fHistDeltaRDDB","Angle D^{*}-D^{0} (Background);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1299  fHistDeltaRDpisB = new TH2F("fHistDeltaRDpisB","Angle D^{*}-#pi_{soft} (Background);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1300  fHistDeltaRDpiB = new TH2F("fHistDeltaRDpiB","Angle D^{*}-#pi (Background);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1301  fHistDeltaRDKB = new TH2F("fHistDeltaRDKB","Angle D^{*}-K (Background);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1302 
1303  fOutput->Add(fHistAlphaDDS);
1304  fOutput->Add(fHistAlphaDpisS);
1305  fOutput->Add(fHistAlphaDpiS);
1306  fOutput->Add(fHistAlphaDKS);
1307  fOutput->Add(fHistAlphaDDB);
1308  fOutput->Add(fHistAlphaDpisB);
1309  fOutput->Add(fHistAlphaDpiB);
1310  fOutput->Add(fHistAlphaDKB);
1311 
1312  fOutput->Add(fHistDeltaRDDS);
1313  fOutput->Add(fHistDeltaRDpisS);
1314  fOutput->Add(fHistDeltaRDpiS);
1315  fOutput->Add(fHistDeltaRDKS);
1316  fOutput->Add(fHistDeltaRDDB);
1317  fOutput->Add(fHistDeltaRDpisB);
1318  fOutput->Add(fHistDeltaRDpiB);
1319  fOutput->Add(fHistDeltaRDKB);
1320  }
1321  else if (fCandidateType == kD0toKpi) {
1322 
1323  fHistAlphaDpiS = new TH2F("fHistAlphaDpiS","Angle D^{0}-#pi (Signal);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1324  fHistAlphaDKS = new TH2F("fHistAlphaDKS","Angle D^{0}-K (Signal);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1325  fHistAlphaDpiR = new TH2F("fHistAlphaDpiR","Angle D^{0}-#pi (Reflections);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1326  fHistAlphaDKR = new TH2F("fHistAlphaDKR","Angle D^{0}-K (Reflections);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1327 
1328  fHistAlphaDpiB = new TH2F("fHistAlphaDpiB","Angle D^{0}-#pi (Background);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1329  fHistAlphaDKB = new TH2F("fHistAlphaDKB","Angle D^{0}-K (Background);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1330 
1331  fHistDeltaRDpiS = new TH2F("fHistDeltaRDpiS","Angle D^{0}-#pi (Signal);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1332  fHistDeltaRDKS = new TH2F("fHistDeltaRDKS","Angle D^{0}-K (Signal);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1333  fHistDeltaRDpiR = new TH2F("fHistDeltaRDpiR","Angle D^{0}-#pi (Reflections);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1334  fHistDeltaRDKR = new TH2F("fHistDeltaRDKR","Angle D^{0}-K (Reflections);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1335 
1336  fHistDeltaRDpiB = new TH2F("fHistDeltaRDpiB","Angle D^{0}-#pi (Background);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1337  fHistDeltaRDKB = new TH2F("fHistDeltaRDKB","Angle D^{0}-K (Background);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1338 
1339  fOutput->Add(fHistAlphaDpiS);
1340  fOutput->Add(fHistAlphaDKS);
1341  fOutput->Add(fHistAlphaDpiR);
1342  fOutput->Add(fHistAlphaDKR);
1343  fOutput->Add(fHistAlphaDpiB);
1344  fOutput->Add(fHistAlphaDKB);
1345 
1346  fOutput->Add(fHistDeltaRDpiS);
1347  fOutput->Add(fHistDeltaRDKS);
1348  fOutput->Add(fHistDeltaRDpiR);
1349  fOutput->Add(fHistDeltaRDKR);
1350  fOutput->Add(fHistDeltaRDpiB);
1351  fOutput->Add(fHistDeltaRDKB);
1352  }
1353  }
1354 
1355  return kTRUE;
1356 }
1357 
1358 //_______________________________________________________________________________
1359 Float_t AliAnalysisTaskSEDmesonsFilterCJ::DeltaR(AliVParticle *p1, AliVParticle *p2) const
1360 {
1361  // Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1362 
1363  if (!p1 || !p2) return -1;
1364 
1365  Double_t phi1 = p1->Phi();
1366  Double_t eta1 = p1->Eta();
1367  Double_t phi2 = p2->Phi();
1368  Double_t eta2 = p2->Eta();
1369 
1370  Double_t dPhi = phi1 - phi2;
1371  if(dPhi < -(TMath::Pi())/2) dPhi = dPhi + TMath::TwoPi();
1372  if(dPhi > (3*(TMath::Pi()))/2) dPhi = dPhi - TMath::TwoPi();
1373 
1374  Double_t dEta = eta1 - eta2;
1375 
1376  Double_t deltaR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
1377 
1378  return deltaR;
1379 }
1380 
1381 //_______________________________________________________________________________
1383 {
1384  //
1385  // Add event tracks to a collection that already contains the D candidates, excluding the daughters of the D candidates
1386  //
1387 
1388  if (!tracks) return;
1389 
1390  TObjArray allDaughters(10);
1391  allDaughters.SetOwner(kFALSE);
1392 
1393  TIter next(coll);
1394  AliEmcalParticle* emcpart = 0;
1395  while ((emcpart = static_cast<AliEmcalParticle*>(next()))) {
1396  AliAODRecoDecay* reco = dynamic_cast<AliAODRecoDecay*>(emcpart->GetTrack());
1397  AliDebug(2, Form("Found a D meson candidtate with pT = %.3f, eta = %.3f, phi = %.3f", reco->Pt(), reco->Eta(), reco->Phi()));
1398  if (reco) AddDaughters(reco, allDaughters);
1399  }
1400 
1401  tracks->ResetCurrentID();
1402  AliVTrack* track = 0;
1403  Int_t n = coll->GetEntriesFast();
1404 
1405  while ((track = static_cast<AliVTrack*>(tracks->GetNextAcceptParticle()))) {
1406  if(fUseMCInfo && fUsePythia){
1407  AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(track);
1408  bool isInj = IsTrackInjected(aodtrack, fMCHeader, fMCarray);
1409  if(!isInj) continue;
1410  }
1411 
1412  if (allDaughters.Remove(track) == 0) {
1413  if(fUseRejTracks){
1414  if(fRan->Rndm() < fTrackIneff) continue;
1415  }
1416 
1417  new ((*coll)[n]) AliEmcalParticle(track);
1418  n++;
1419  AliDebug(2, Form("Track %d (pT = %.3f, eta = %.3f, phi = %.3f) is included", tracks->GetCurrentID(), track->Pt(), track->Eta(), track->Phi()));
1420  }
1421  else {
1422  AliDebug(2, Form("Track %d (pT = %.3f, eta = %.3f, phi = %.3f) is excluded", tracks->GetCurrentID(), track->Pt(), track->Eta(), track->Phi()));
1423  }
1424  }
1425 }
1426 
1427 //_______________________________________________________________________________
1429 {
1430  // Add all the dauthers of cand in an array. Follows all the decay cascades.
1431 
1432  Int_t n = cand->GetNDaughters();
1433 
1434  //Printf("AddDaughters: the number of dauhters is %d", n);
1435 
1436  Int_t ntot = 0;
1437  Double_t pt = 0;
1438  for (Int_t i = 0; i < n; i++) {
1439  AliVTrack* track = dynamic_cast<AliVTrack*>(cand->GetDaughter(i));
1440  if (!track) continue;
1441 
1442  AliAODRecoDecay* cand2 = dynamic_cast<AliAODRecoDecay*>(track);
1443 
1444  if (cand2) {
1445  //Printf("Daughter pT = %.3f --> ", track->Pt());
1446  pt += AddDaughters(cand2, daughters);
1447  }
1448  else {
1449  if (!track->InheritsFrom("AliAODTrack")) {
1450  Printf("Warning: One of the daughters is not of type 'AliAODTrack' nor 'AliAODRecoDecay'.");
1451  continue;
1452  }
1453  //Printf("Daughter pT = %.3f", track->Pt());
1454  daughters.AddLast(track);
1455  pt += track->Pt();
1456  ntot++;
1457  }
1458  }
1459 
1460  //Printf("Total pt of the daughters = %.3f", pt);
1461 
1462  return pt;
1463 }
1464 //_______________________________________________________________________________
1466 {
1467  //
1468  // Add event tracks to a collection that already contains the D candidates, excluding the daughters of the D candidates
1469  //
1470 
1471  if (!mctracks) return;
1472 
1473  TObjArray allMCDaughters(10);
1474  allMCDaughters.SetOwner(kFALSE);
1475 
1476  AliAODMCParticle* mcD = (AliAODMCParticle*)coll->At(0);
1477 
1478  AddMCDaughters(mcD,allMCDaughters,fMCarray);
1479 
1480  mctracks->ResetCurrentID();
1481  AliAODMCParticle* mcpart = 0;
1482  Int_t n = coll->GetEntriesFast();
1483  while ((mcpart = static_cast<AliAODMCParticle*>(mctracks->GetNextAcceptParticle()))) {
1484  if(TMath::Abs(mcpart->Charge())==0) continue;
1485  if(fUsePythia){
1486  bool isInj = IsMCTrackInjected(mcpart, fMCHeader, fMCarray);
1487  if(!isInj) continue;
1488  }
1489 
1490 
1491  if (allMCDaughters.Remove(mcpart) == 0) {
1492  new ((*coll)[n]) AliAODMCParticle(*mcpart);
1493  n++;
1494  AliDebug(2, Form("Track %d (pT = %.3f, eta = %.3f, phi = %.3f) is included", mctracks->GetCurrentID(), mcpart->Pt(), mcpart->Eta(), mcpart->Phi()));
1495  }
1496  else {
1497  AliDebug(2, Form("Track %d (pT = %.3f, eta = %.3f, phi = %.3f) is excluded", mctracks->GetCurrentID(), mcpart->Pt(), mcpart->Eta(), mcpart->Phi()));
1498  }
1499  }
1500 }
1501 
1502 //_______________________________________________________________________________
1503 Double_t AliAnalysisTaskSEDmesonsFilterCJ::AddMCDaughters(AliAODMCParticle* mcDmeson, TObjArray& mcdaughters, TClonesArray* mcArray)
1504 {
1505  // Add all the dauthers of cand in an array. Follows all the decay cascades.
1506 
1507  Int_t n = mcDmeson->GetNDaughters();
1508  Int_t nD0 = mcDmeson->GetDaughter(0); // get label of the first daughter
1509  //Printf("AddDaughters: the number of dauhters is %d", n);
1510  Double_t pt = 0;
1511 
1512  for (Int_t i = 0; i < n; i++) {
1513  AliAODMCParticle* DDaughter = static_cast<AliAODMCParticle*>(mcArray->At(nD0+i));
1514  if (!DDaughter) continue;
1515 
1516  if (DDaughter->GetNDaughters()>0) {
1517  //Printf("Daughter pT = %.3f --> ", track->Pt());
1518  pt += AddMCDaughters(DDaughter, mcdaughters, mcArray);
1519  mcdaughters.AddLast(DDaughter);
1520  }
1521  else {
1522  //Printf("Daughter pT = %.3f", track->Pt());
1523  mcdaughters.AddLast(DDaughter);
1524  pt += DDaughter->Pt();
1525  }
1526  }
1527 
1528  //Printf("Total pt of the daughters = %.3f", pt);
1529 
1530  return pt;
1531 }
1532 //_________________________________________________________________________________________________
1533 Int_t AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin(AliAODRecoDecay* cand, TClonesArray* mcArray)
1534 {
1535  // Checks whether the mother of the D meson candidate comes from a charm or a bottom quark.
1536 
1537  if (!mcArray) return -1;
1538 
1539  Int_t labDau0 = static_cast<AliVTrack*>(cand->GetDaughter(0))->GetLabel();
1540  if (labDau0 < 0) return -1;
1541 
1542  AliAODMCParticle* part = static_cast<AliAODMCParticle*>(mcArray->At(labDau0));
1543  return CheckOrigin(part, mcArray);
1544 }
1545 
1546 //_________________________________________________________________________________________________
1548 {
1549  // Checks whether the mother of the D meson candidate comes from a charm or a bottom quark.
1550 
1551  if (!stack) return -1;
1552 
1553  Int_t labDau0 = static_cast<AliVTrack*>(cand->GetDaughter(0))->GetLabel();
1554  if (labDau0 < 0) return -1;
1555 
1556  return CheckOrigin(labDau0, stack);
1557 }
1558 
1559 //_________________________________________________________________________________________________
1560 Int_t AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin(AliAODMCParticle* part, TClonesArray* mcArray)
1561 {
1562  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1563 
1564  if (!part) return -1;
1565  if (!mcArray) return -1;
1566 
1567  Int_t pdgGranma = 0;
1568  Int_t mother = part->GetMother();
1569  Int_t istep = 0;
1570  Int_t abspdgGranma = 0;
1571  Bool_t isFromB = kFALSE;
1572  Bool_t isQuarkFound = kFALSE;
1573 
1574  while (mother >= 0) {
1575  istep++;
1576  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1577  if (mcGranma != 0) {
1578  pdgGranma = mcGranma->GetPdgCode();
1579  abspdgGranma = TMath::Abs(pdgGranma);
1580  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1581  isFromB = kTRUE;
1582  }
1583 
1584  if (abspdgGranma == 4 || abspdgGranma == 5) isQuarkFound = kTRUE;
1585  mother = mcGranma->GetMother();
1586  }
1587  else {
1588  ::Error("AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1589  break;
1590  }
1591  }
1592 
1593  if (isQuarkFound) {
1594  if (isFromB) {
1595  return kFromBottom;
1596  }
1597  else {
1598  return kFromCharm;
1599  }
1600  }
1601  else {
1602  return kQuarkNotFound;
1603  }
1604 }
1605 
1606 //_________________________________________________________________________________________________
1608 {
1609  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1610 
1611  if (!stack) return -1;
1612 
1613  TParticle* part = stack->Particle(ipart);
1614  if (!part) return -1;
1615 
1616  Int_t pdgGranma = 0;
1617  Int_t mother = part->GetFirstMother();
1618  Int_t istep = 0;
1619  Int_t abspdgGranma = 0;
1620  Bool_t isFromB = kFALSE;
1621  Bool_t isQuarkFound = kFALSE;
1622 
1623  while (mother >= 0) {
1624  istep++;
1625  TParticle* mcGranma = stack->Particle(mother);
1626  if (mcGranma != 0) {
1627  pdgGranma = mcGranma->GetPdgCode();
1628  abspdgGranma = TMath::Abs(pdgGranma);
1629  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1630  isFromB = kTRUE;
1631  }
1632 
1633  if (abspdgGranma == 4 || abspdgGranma == 5) isQuarkFound = kTRUE;
1634  mother = mcGranma->GetFirstMother();
1635  }
1636  else {
1637  ::Error("AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1638  break;
1639  }
1640  }
1641 
1642  if (isQuarkFound) {
1643  if (isFromB) {
1644  return kFromBottom;
1645  }
1646  else {
1647  return kFromCharm;
1648  }
1649  }
1650  else {
1651  return kQuarkNotFound;
1652  }
1653 }
1654 
1655 //_________________________________________________________________________________________________
1656 Int_t AliAnalysisTaskSEDmesonsFilterCJ::CheckDecayChannel(AliAODMCParticle* part, TClonesArray* mcArray)
1657 {
1658  // Determine the decay channel
1659 
1660  if (!part) return -1;
1661  if (!mcArray) return -1;
1662 
1664 
1665  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1666 
1667  if (part->GetNDaughters() == 2) {
1668 
1669  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1670  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1671 
1672  if (!d1 || !d2) {
1673  return decay;
1674  }
1675 
1676  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1677  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1678 
1679  if (absPdgPart == 421) { // D0 -> K pi
1680 
1681  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1682  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1683  decay = kDecayD0toKpi;
1684  }
1685  }
1686 
1687  if (absPdgPart == 413) { // D* -> D0 pi
1688 
1689  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1690  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1691  if (D0decay == kDecayD0toKpi) {
1692  decay = kDecayDStartoKpipi;
1693  }
1694  }
1695 
1696  if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1697  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1698  if (D0decay == kDecayD0toKpi) {
1699  decay = kDecayDStartoKpipi;
1700  }
1701  }
1702  }
1703  }
1704 
1705  return decay;
1706 }
1707 
1708 //_________________________________________________________________________________________________
1710 {
1711  // Determine the decay channel
1712 
1713  if (!stack) return -1;
1714 
1715  TParticle* part = stack->Particle(ipart);
1716 
1717  if (!part) return -1;
1718 
1720 
1721  if (part->GetNDaughters() == 2) {
1722 
1723  Int_t id1 = part->GetDaughter(0);
1724  Int_t id2 = part->GetDaughter(1);
1725 
1726  TParticle* d1 = stack->Particle(id1);
1727  TParticle* d2 = stack->Particle(id2);
1728 
1729  if (!d1 || !d2) {
1730  return decay;
1731  }
1732 
1733  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1734  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1735 
1736 
1737  if (part->GetPdgCode() == 421) { // D0 -> K pi
1738 
1739  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1740  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1741  decay = kDecayD0toKpi;
1742  }
1743  }
1744 
1745  if (part->GetPdgCode() == 413) { // D* -> D0 pi
1746 
1747  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1748  Int_t D0decay = CheckDecayChannel(id1, stack);
1749  if (D0decay == kDecayD0toKpi) {
1750  decay = kDecayDStartoKpipi;
1751  }
1752  }
1753 
1754  if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1755  Int_t D0decay = CheckDecayChannel(id2, stack);
1756  if (D0decay == kDecayD0toKpi) {
1757  decay = kDecayDStartoKpipi;
1758  }
1759  }
1760  }
1761  }
1762 
1763  return decay;
1764 
1765 }
1766 
1767 
1768 //_____________________________________________________________________
1769 void AliAnalysisTaskSEDmesonsFilterCJ::GetTrackPrimaryGenerator(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){
1770 
1772 
1773  Int_t lab=TMath::Abs(track->GetLabel());
1774  nameGen=AliVertexingHFUtils::GetGenerator(lab,header);
1775 
1776  // Int_t countControl=0;
1777 
1778  while(nameGen.IsWhitespace()){
1779  AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);
1780  if(!mcpart){
1781  printf("AliAnalysisTaskMultCheck::IsTrackInjected - BREAK: No valid AliAODMCParticle at label %i\n",lab);
1782  break;
1783  }
1784  Int_t mother = mcpart->GetMother();
1785  if(mother<0){
1786  printf("AliAnalysisTaskMultCheck::IsTrackInjected - BREAK: Reached primary particle without valid mother\n");
1787  break;
1788  }
1789  lab=mother;
1790  nameGen=AliVertexingHFUtils::GetGenerator(mother,header);
1791  // countControl++;
1792  // if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops
1793  // printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Protection from infinite loop active\n");
1794  // break;
1795  // }
1796  }
1797 
1798  return;
1799 }
1800 //----------------------------------------------------------------------
1801 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::IsTrackInjected(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC){
1803  TString nameGen;
1804  GetTrackPrimaryGenerator(track,header,arrayMC,nameGen);
1805 
1806  if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
1807 
1808  return kTRUE;
1809 }
1810 
1811 
1812 //_____________________________________________________________________
1813 void AliAnalysisTaskSEDmesonsFilterCJ::GetMCTrackPrimaryGenerator(AliAODMCParticle *track,AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){
1814 
1816 
1817  Int_t lab=TMath::Abs(track->GetLabel());
1818  nameGen=AliVertexingHFUtils::GetGenerator(lab,header);
1819 
1820  // Int_t countControl=0;
1821 
1822  while(nameGen.IsWhitespace()){
1823  AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);
1824  if(!mcpart){
1825  printf("AliAnalysisTaskMultCheck::IsTrackInjected - BREAK: No valid AliAODMCParticle at label %i\n",lab);
1826  break;
1827  }
1828  Int_t mother = mcpart->GetMother();
1829  if(mother<0){
1830  printf("AliAnalysisTaskMultCheck::IsTrackInjected - BREAK: Reached primary particle without valid mother\n");
1831  break;
1832  }
1833  lab=mother;
1834  nameGen=AliVertexingHFUtils::GetGenerator(mother,header);
1835  // countControl++;
1836  // if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops
1837  // printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Protection from infinite loop active\n");
1838  // break;
1839  // }
1840  }
1841 
1842  return;
1843 }
1844 //----------------------------------------------------------------------
1845 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::IsMCTrackInjected(AliAODMCParticle *track,AliAODMCHeader *header,TClonesArray *arrayMC){
1847  TString nameGen;
1848  GetMCTrackPrimaryGenerator(track,header,arrayMC,nameGen);
1849 
1850  if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
1851 
1852  return kTRUE;
1853 }
void ProcessDstar(AliAODRecoCascadeHF *dstar, Int_t isSelected)
Int_t pdg
static Int_t CheckOrigin(AliAODRecoDecay *cand, TClonesArray *mcArray)
virtual AliVParticle * GetNextAcceptParticle()
double Double_t
Definition: External.C:58
Definition: External.C:236
TClonesArray * fCombinedDmesonsBkg
contains candidates selected by AliRDHFCuts and the rest of the event tracks
Double_t DeltaInvMass() const
void ProcessD0(AliAODRecoDecayHF2Prong *charmCand, Int_t isSelected)
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
Base task in the EMCAL framework.
TClonesArray * fArrayDStartoD0pi
Random number generator.
Double_t mass
static Int_t CheckMatchingAODdeltaAODevents()
static TString GetGenerator(Int_t label, AliAODMCHeader *header)
void FillD0MCTruthKinHistos(AliAODRecoDecayHF2Prong *charmCand, Int_t isSelected, Int_t isD0)
static Double_t AddDaughters(AliAODRecoDecay *cand, TObjArray &daughters)
void AddEventTracks(TClonesArray *coll, AliParticleContainer *tracks)
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
AliVTrack * GetTrack() const
AliStack * stack
void GetMCTrackPrimaryGenerator(AliAODMCParticle *track, AliAODMCHeader *header, TClonesArray *arrayMC, TString &nameGen)
void FillDStarMCTruthKinHistos(AliAODRecoCascadeHF *dstar, Int_t isSelected, Int_t isDstar)
static Double_t AddMCDaughters(AliAODMCParticle *mcDmeson, TObjArray &mcdaughters, TClonesArray *mcArray)
Container for particles within the EMCAL framework.
Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Int_t fNCand
contains MC D0 and MC event particles
Bool_t IsTrackInjected(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC)
void FillDstarSideBands(AliAODRecoCascadeHF *dstar)
TClonesArray * fSideBandArray
contains candidates selected by AliRDHFCuts
int Int_t
Definition: External.C:63
Definition: External.C:204
static Int_t CheckDecayChannel(AliAODMCParticle *part, TClonesArray *mcArray)
float Float_t
Definition: External.C:68
const Double_t ptmax
AliAODTrack * GetBachelor() const
TClonesArray * fMCCombinedDmesons
contains bkg candidates selected by AliRDHFCuts and the rest of the event tracks
const Double_t ptmin
Bool_t SetD0WidthForDStar(Int_t nptbins, Float_t *width)
void ExecOnce()
Perform steps needed to initialize the analysis.
Bool_t IsEventSelected(AliVEvent *event)
decay
Definition: HFPtSpectrum.C:41
AliEmcalList * fOutput
!output list
Int_t fNSBCand
number of selected D candidates already added to fCandidateArray
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:292
virtual void ExecOnce()
Perform steps needed to initialize the analysis.
Bool_t IsMCTrackInjected(AliAODMCParticle *track, AliAODMCHeader *header, TClonesArray *arrayMC)
Double_t InvMassD0() const
TH1 * fHistStat
number of selected side-band D candidates already added to fSideBandArray
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:249
Float_t DeltaR(AliVParticle *p1, AliVParticle *p2) const
void AddObjectToEvent(TObject *obj, Bool_t attempt=kFALSE)
Add object to event.
void UserCreateOutputObjects()
Main initialization function on the worker.
bool Bool_t
Definition: External.C:53
AliAODRecoDecayHF2Prong * Get2Prong() const
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:312
Int_t PtBin(Double_t pt) const
void AddMCEventTracks(TClonesArray *coll, AliParticleContainer *mctracks)
Bool_t fNeedEmcalGeom
whether or not the task needs the emcal geometry
Int_t nptbins
void GetTrackPrimaryGenerator(AliAODTrack *track, AliAODMCHeader *header, TClonesArray *arrayMC, TString &nameGen)
TClonesArray * fCombinedDmesons
contains candidates selected by AliRDHFCuts::IsSelected(kTracks), to be used for side bands (DStar ca...