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