AliPhysics  857879d (857879d)
 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  fHistStat->Fill(3);
549 
550  if (fNCand==fAnalyseCand){
551 
552  new ((*fMCCombinedDmesons)[0]) AliAODMCParticle(*charmPart);
553 
554  AliAODRecoDecayHF2Prong* charmCand = 0;
555  AliAODRecoCascadeHF* dstar = 0;
556 
557  // loop over reco D candidates to find a match to MC
558  for (Int_t icharm = 0; icharm < nD; icharm++) {
559 
560  charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fArrayDStartoD0pi->At(icharm)); // D candidates
561  if (!charmCand) continue;
562  if(!(vHF->FillRecoCand(fAodEvent,charmCand))) continue;
563 
564  Int_t nprongs = charmCand->GetNProngs();
565  AliDebug(2, Form("Candidate is %d, and nprongs = %d", fCandidateType, nprongs));
566 
567  if (fCandidateType == kDstartoKpipi) {
568  dstar = dynamic_cast<AliAODRecoCascadeHF*>(charmCand);
569  if (!dstar) {
570  Error("AliAnalysisTaskSEDmesonsFilterCJ::UserExec","Candidate type is D* but object type is wrong (should be AliAODRecoCascadeHF)");
571  continue;
572  }
573  }
574 
575  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
576  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
577 
578  Int_t mcLabel = NULL;
579  if (fCandidateType == kDstartoKpipi) mcLabel = dstar->MatchToMC(413, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCarray);
580  else mcLabel = charmCand->MatchToMC(421, fMCarray, fNProngs, fPDGdaughters);
581 
582  if(mcLabel == iMCcharm) break;
583  }
584 
585  if (!charmCand) break;
586  if (fCandidateType == kDstartoKpipi && !dstar) break;
587 
588  fHistStat->Fill(4);
589 
590  // region of interest + cuts
591  if (!fCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(pdgMeson))) break;
592 
593  //candidate selected by cuts and PID
594  Int_t isSelected = 0;
595  isSelected = fCuts->IsSelected(charmCand, AliRDHFCuts::kAll, fAodEvent); //selected
596  if (!isSelected) break;
597 
598  fHistStat->Fill(5);
599 
600  if (fCandidateType == kDstartoKpipi) {
601  AliAODRecoDecayHF2Prong* D0fromDstar = dstar->Get2Prong();
602  fHistInvMassS->Fill(dstar->DeltaInvMass());
603  fHistImpParS->Fill(dstar->Getd0Prong(0), dstar->PtProng(0)); //bachelor
604  fHistImpParS->Fill(D0fromDstar->Getd0Prong(0), D0fromDstar->PtProng(0));
605  fHistImpParS->Fill(D0fromDstar->Getd0Prong(1), D0fromDstar->PtProng(1));
606  }
607  else {
608  fHistImpParS->Fill(charmCand->Getd0Prong(0), charmCand->PtProng(0));
609  fHistImpParS->Fill(charmCand->Getd0Prong(1), charmCand->PtProng(1));
610  fHistInvMassS->Fill(charmCand->InvMassD0());
611  fHistInvMassS->Fill(charmCand->InvMassD0bar());
612  }
613 
614  if (fCandidateType == kDstartoKpipi) {
615  Int_t isDstar = charmPart == NULL ? 0 : 1;
616  //fill histograms of kinematics, using MC truth
617  FillDStarMCTruthKinHistos(dstar, isSelected, isDstar);
618 
619  new ((*fCandidateArray)[0]) AliAODRecoCascadeHF(*dstar);
620  new ((*fCombinedDmesons)[0]) AliEmcalParticle(dstar);
621 
622  }
623  else {
624 
625  Int_t isD0 = 0;
626  Int_t pdgCode = charmPart->GetPdgCode();
627  if (pdgCode == 421) { isD0 = 1; }
628  else if (pdgCode == -421) { isD0 = -1; }
629  //fill histograms of kinematics, using MC truth
630  FillD0MCTruthKinHistos(charmCand, isSelected, isD0);
631 
632  new ((*fCandidateArray)[0]) AliAODRecoDecayHF2Prong(*charmCand);
633  new ((*fCombinedDmesons)[0]) AliEmcalParticle(charmCand);
634  }
635 
636  break;
637  }
638  else { fNCand++; }
639 
640 }
641 }
642 // for data, or MC without RM or efficiency studies
643 else {
644  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
645  Int_t isSelected = 0;
646 
647  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fArrayDStartoD0pi->At(icharm)); // D candidates
648  if (!charmCand) continue;
649  if(!(vHF->FillRecoCand(fAodEvent,charmCand))) continue;
650 
651  Int_t nprongs = charmCand->GetNProngs();
652  AliDebug(2, Form("Candidate is %d, and nprongs = %d", fCandidateType, nprongs));
653 
654  AliAODRecoCascadeHF* dstar = 0;
655 
656  if (fCandidateType == kDstartoKpipi) {
657  dstar = dynamic_cast<AliAODRecoCascadeHF*>(charmCand);
658  if (!dstar) {
659  Error("AliAnalysisTaskSEDmesonsFilterCJ::UserExec","Candidate type is D* but object type is wrong (should be AliAODRecoCascadeHF)");
660  continue;
661  }
662  }
663 
664  // region of interest + cuts
665  if (!fCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(pdgMeson))) continue;
666 
668  FillDstarSideBands(dstar);
669  }
670 
671 
672  //candidate selected by cuts and PID
673  isSelected = fCuts->IsSelected(charmCand, AliRDHFCuts::kAll, fAodEvent); //selected
674  if (!isSelected) continue;
675 
676  if (fCandidateType == kDstartoKpipi) {
677  ProcessDstar(dstar, isSelected);
678  }
679  else {
680  ProcessD0(charmCand, isSelected);
681  }
682  } // end of D cand loop
683 }
684 
685 
686  delete vHF;
687 
688 
689  AliDebug(2, "Loop done");
690 
691  if (fCombineDmesons) {
692  if (fCombinedDmesons->GetEntriesFast() > 0) {
694  }
695 
696  if (fMCCombinedDmesons->GetEntriesFast() > 0) {
698  }
699 
700  if (fCombinedDmesonsBkg->GetEntriesFast() > 0) {
702  }
703  }
704 
705  fHistNCandEv->Fill(fCandidateArray->GetEntriesFast());
707  Int_t nsbcand = fSideBandArray->GetEntriesFast();
708  fHistStat->Fill(4, nsbcand);
709  fHistNSBCandEv->Fill(nsbcand);
710  }
711 
712  PostData(1, fOutput);
713  PostData(3, fCandidateArray);
714  PostData(4, fSideBandArray);
715  PostData(5, fCombinedDmesons);
716  PostData(6, fCombinedDmesonsBkg);
717  PostData(7, fMCCombinedDmesons);
718 
719  AliDebug(2, "Exiting method");
720 
721  return kTRUE;
722 }
723 
724 //_______________________________________________________________________________
726 {
727  // Process the D0 candidate.
728 
729  // For MC analysis look for the AOD MC particle associated with the charm candidate
730  AliAODMCParticle* charmPart = 0;
731  if (fUseMCInfo) {
732  Int_t mcLabel = charmCand->MatchToMC(421, fMCarray, fNProngs, fPDGdaughters);
733 
734  if (mcLabel >= 0) {
735  charmPart = static_cast<AliAODMCParticle*>(fMCarray->At(mcLabel));
736  }
737  }
738 
739  if (charmPart) {
740 
741  Int_t origin = CheckOrigin(charmPart, fMCarray);
742  if (origin < 0) return;
743 
744  if (fRejectQuarkNotFound && origin == kQuarkNotFound) {
745  fHistStat->Fill(6);
746  return;
747  }
748  if (fRejectDfromB && origin == kFromBottom) {
749  fHistStat->Fill(7);
750  return;
751  }
752  if (fKeepOnlyDfromB && origin != kFromBottom) {
753  fHistStat->Fill(8);
754  return;
755  }
756 
757  fHistStat->Fill(2);
758  }
759  else {
760  fHistStat->Fill(5);
761  }
762 
763  // For MC background fill fSideBandArray
764  if (fUseMCInfo && !charmPart) {
765  if (fUseReco) {
766  if(!fMultCand) new ((*fSideBandArray)[fNSBCand]) AliAODRecoDecayHF2Prong(*charmCand);
767  else if(fNSBCand==fAnalyseCand) new ((*fSideBandArray)[0]) AliAODRecoDecayHF2Prong(*charmCand);
768 
769  if (fCombineDmesons) {
770  if(!fMultCand) new ((*fCombinedDmesonsBkg)[fNSBCand]) AliEmcalParticle(charmCand);
771  else if(fNSBCand==fAnalyseCand) new ((*fCombinedDmesonsBkg)[0]) AliEmcalParticle(charmCand);
772  }
773 
774  fHistImpParB->Fill(charmCand->Getd0Prong(0), charmCand->PtProng(0));
775  fHistImpParB->Fill(charmCand->Getd0Prong(1), charmCand->PtProng(1));
776  fHistInvMassB->Fill(charmCand->InvMassD0());
777  fHistInvMassB->Fill(charmCand->InvMassD0bar());
778 
779  fNSBCand++;
780  }
781  }
782  // For data or MC signal fill fCandidateArray
783  else {
784  // For data or MC with the requirement fUseReco fill with candidates
785  if (fUseReco) {
786  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODRecoDecayHF2Prong(*charmCand);
787  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODRecoDecayHF2Prong(*charmCand);
788 
789  if (fCombineDmesons) {
790  if(!fMultCand) new ((*fCombinedDmesons)[fNCand]) AliEmcalParticle(charmCand);
791  else if(fNCand==fAnalyseCand)
792  {
793  new ((*fCombinedDmesons)[0]) AliEmcalParticle(charmCand);
794  if(charmPart) new ((*fMCCombinedDmesons)[0]) AliAODMCParticle(*charmPart);
795  }
796  }
797 
798  fHistImpParS->Fill(charmCand->Getd0Prong(0), charmCand->PtProng(0));
799  fHistImpParS->Fill(charmCand->Getd0Prong(1), charmCand->PtProng(1));
800  fHistInvMassS->Fill(charmCand->InvMassD0());
801  fHistInvMassS->Fill(charmCand->InvMassD0bar());
802  }
803  // For MC with requirement particle level fill with AliAODMCParticle
804  else {
805  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODMCParticle(*charmPart);
806  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODMCParticle(*charmPart);
807  }
808  fHistStat->Fill(3);
809  fNCand++;
810  }
811 
812  // Now filling some more histograms
813 
814  // mass vs pt
815  if (isSelected == 1 || isSelected == 3) fHistInvMassPtD->Fill(charmCand->InvMassD0(), charmCand->Pt());
816  if (isSelected >= 2) fHistInvMassPtD->Fill(charmCand->InvMassD0bar(), charmCand->Pt());
817 
818  if (fUseMCInfo) { //fill histograms of kinematics, using MC truth
819  Int_t isD0 = 0;
820  if (charmPart) {
821  Int_t pdgCode = charmPart->GetPdgCode();
822 
823  if (pdgCode == 421) {
824  isD0 = 1;
825  }
826  else if (pdgCode == -421) {
827  isD0 = -1;
828  }
829  else {
830  AliDebug(2, "Not a D0/D0bar meson!");
831  return;
832  }
833  }
834 
835  FillD0MCTruthKinHistos(charmCand, isSelected, isD0);
836  }
837 }
838 
839 //_______________________________________________________________________________
841 {
842  // Process the D* candidate.
843 
844  AliDebug(2, "Entering method");
845 
846  // For MC analysis look for the AOD MC particle associated with the charm candidate
847  AliAODMCParticle* charmPart = 0;
848  if (fUseMCInfo) {
849  //D* and D0 prongs needed to MatchToMC method
850  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
851  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
852 
853  Int_t mcLabel = dstar->MatchToMC(413, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCarray);
854 
855  if (mcLabel >= 0) {
856  charmPart = static_cast<AliAODMCParticle*>(fMCarray->At(mcLabel));
857  }
858  }
859 
860  if (charmPart) {
861 
862  Int_t origin = CheckOrigin(charmPart, fMCarray);
863  if (origin < 0) return;
864 
865  if (fRejectQuarkNotFound && origin == kQuarkNotFound) {
866  fHistStat->Fill(6);
867  return;
868  }
869  if (fRejectDfromB && origin == kFromBottom) {
870  fHistStat->Fill(7);
871  return;
872  }
873  if (fKeepOnlyDfromB && origin != kFromBottom) {
874  fHistStat->Fill(8);
875  return;
876  }
877 
878  fHistStat->Fill(2);
879  }
880  else {
881  fHistStat->Fill(5);
882  }
883 
884  AliAODRecoDecayHF2Prong* D0fromDstar = dstar->Get2Prong();
885 
886  // For MC background fill fSideBandArray
887  if (fUseMCInfo && !charmPart) {
888  if (fUseReco) {
889  if(!fMultCand) new ((*fSideBandArray)[fNSBCand]) AliAODRecoCascadeHF(*dstar);
890  else if(fNSBCand==fAnalyseCand) new ((*fSideBandArray)[0]) AliAODRecoCascadeHF(*dstar);
891 
892  if (fCombineDmesons) {
893  if(!fMultCand) new ((*fCombinedDmesonsBkg)[fNSBCand]) AliEmcalParticle(dstar);
894  else if(fNSBCand==fAnalyseCand) new ((*fCombinedDmesonsBkg)[0]) AliEmcalParticle(dstar);
895  }
896 
897  fHistInvMassB->Fill(dstar->DeltaInvMass());
898  fHistImpParB->Fill(dstar->Getd0Prong(0), dstar->PtProng(0)); //bachelor
899  fHistImpParB->Fill(D0fromDstar->Getd0Prong(0), D0fromDstar->PtProng(0));
900  fHistImpParB->Fill(D0fromDstar->Getd0Prong(1), D0fromDstar->PtProng(1));
901 
902  fNSBCand++;
903  }
904  }
905  // For data and MC signal fill fCandidateArray
906  else {
907  // For data or MC signal with the requirement fUseReco fill with candidates
908  if (fUseReco) {
909  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODRecoCascadeHF(*dstar);
910  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODRecoCascadeHF(*dstar);
911 
912  if (fCombineDmesons) {
913  if(!fMultCand) new ((*fCombinedDmesons)[fNCand]) AliEmcalParticle(dstar);
914  else if(fNCand==fAnalyseCand)
915  {
916  new ((*fCombinedDmesons)[0]) AliEmcalParticle(dstar);
917  if(charmPart) new ((*fMCCombinedDmesons)[0]) AliAODMCParticle(*charmPart);
918  }
919  }
920  }
921  // For MC signal with requirement particle level fill with AliAODMCParticle
922  else {
923  if(!fMultCand) new ((*fCandidateArray)[fNCand]) AliAODMCParticle(*charmPart);
924  else if(fNCand==fAnalyseCand) new ((*fCandidateArray)[0]) AliAODMCParticle(*charmPart);
925  }
926  fNCand++;
927 
928  fHistInvMassS->Fill(dstar->DeltaInvMass());
929  fHistImpParS->Fill(dstar->Getd0Prong(0), dstar->PtProng(0)); //bachelor
930  fHistImpParS->Fill(D0fromDstar->Getd0Prong(0), D0fromDstar->PtProng(0));
931  fHistImpParS->Fill(D0fromDstar->Getd0Prong(1), D0fromDstar->PtProng(1));
932  fHistStat->Fill(3);
933  }
934 
935  // Now filling some more histograms.
936 
937  // select D* in the D0 window.
938  // In the cut object window is loose to allow for side bands
939 
940  // retrieve the corresponding pt bin for the candidate
941  // and set the expected D0 width (x3)
942 
943  Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
944  Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
945  Double_t invMassPDG = mPDGDstar - mPDGD0;
946 
947  Double_t ptD = dstar->Pt();
948  Int_t ptDbin = fCuts->PtBin(ptD);
949  if (ptDbin < 0 || ptDbin >= fCuts->GetNPtBins()) {
950  AliError(Form("Pt %.3f out of bounds", ptD));
951  return;
952  }
953 
954  AliDebug(1, Form("Pt bin %d and sigma D0 %.4f", ptDbin, fSigmaD0[ptDbin]));
955  //consider the D* candidates only if the mass of the D0 is in 3 sigma wrt the PDG value
956  if ((dstar->InvMassD0()>=(mPDGD0-3.*fSigmaD0[ptDbin])) && (dstar->InvMassD0()<=(mPDGD0+3.*fSigmaD0[ptDbin]))) {
957  AliDebug(2, "D0 mass within 3 sigma");
958 
959  // D* delta mass
960  AliDebug(2, "Filling invariant mass vs pt histogram");
961  fHistInvMassPtD->Fill(dstar->DeltaInvMass(), ptD); // 2 D slice for pt bins
962 
963  // Soft pion pt for good candidates
964  Double_t invmassDelta = dstar->DeltaInvMass();
965  if (TMath::Abs(invmassDelta - invMassPDG) < 0.0021) {
966  AliDebug(2, "Filling pion pt histogram");
967  //softpion from D* decay
968  AliAODTrack *softPionTrack = static_cast<AliAODTrack*>(dstar->GetBachelor());
969  if (softPionTrack) fHistPtPion->Fill(softPionTrack->Pt());
970  }
971  }
972 
973  if (fUseMCInfo) { //fill histograms of kinematics, using MC truth
974  Int_t isDstar = charmPart == NULL ? 0 : 1;
975  FillDStarMCTruthKinHistos(dstar, isDstar, isSelected);
976  }
977 
978  AliDebug(2, "Exiting method");
979 }
980 
981 //_______________________________________________________________________________
983 {
984  // Fill some histogram with kinematic information of the D0 candidate.
985 
986  Double_t ptD = charmCand->Pt();
987  Double_t aD = charmCand->Phi();
988  Double_t adaugh[2] = {charmCand->PhiProng(0), charmCand->PhiProng(1)};
989  AliAODTrack* p0 = static_cast<AliAODTrack*>(charmCand->GetDaughter(0));
990  AliAODTrack* p1 = static_cast<AliAODTrack*>(charmCand->GetDaughter(1));
991  Float_t dR0 = DeltaR(charmCand, p0);
992  Float_t dR1 = DeltaR(charmCand, p1);
993 
994  if (isD0 == 0) { //background
995  if (isSelected == 1 || isSelected == 3) { // selected as D0
996  fHistAlphaDKB->Fill(aD-adaugh[0],ptD);
997  fHistAlphaDpiB->Fill(aD-adaugh[1],ptD);
998 
999  fHistDeltaRDKB->Fill(dR0,ptD);
1000  fHistDeltaRDpiB->Fill(dR1,ptD);
1001  }
1002  if (isSelected >= 2) { //selected as D0bar
1003  fHistAlphaDpiB->Fill(aD-adaugh[0],ptD);
1004  fHistAlphaDKB->Fill(aD-adaugh[1],ptD);
1005 
1006  fHistDeltaRDpiB->Fill(dR0,ptD);
1007  fHistDeltaRDKB->Fill(dR1,ptD);
1008  }
1009  }
1010  else if (isD0 == 1) { //D0
1011  fHistAlphaDKS->Fill(aD-adaugh[0],ptD);
1012  fHistAlphaDpiS->Fill(aD-adaugh[1],ptD);
1013 
1014  fHistDeltaRDKS->Fill(dR0,ptD);
1015  fHistDeltaRDpiS->Fill(dR1,ptD);
1016 
1017  if (isSelected == 3) { // selected as both D0bar/D0
1018  fHistAlphaDpiR->Fill(aD-adaugh[0],ptD);
1019  fHistAlphaDKR->Fill(aD-adaugh[1],ptD);
1020 
1021  fHistDeltaRDpiR->Fill(dR0,ptD);
1022  fHistDeltaRDKR->Fill(dR1,ptD);
1023  }
1024  }
1025  else if (isD0 == -1) { //D0bar
1026  fHistAlphaDKS->Fill(aD-adaugh[1],ptD);
1027  fHistAlphaDpiS->Fill(aD-adaugh[0],ptD);
1028 
1029  fHistDeltaRDKS->Fill(dR1,ptD);
1030  fHistDeltaRDpiS->Fill(dR0,ptD);
1031 
1032  if (isSelected == 3) { // selected as both D0bar/D0
1033  fHistAlphaDpiR->Fill(aD-adaugh[1],ptD);
1034  fHistAlphaDKR->Fill(aD-adaugh[0],ptD);
1035 
1036  fHistDeltaRDpiR->Fill(dR1, ptD);
1037  fHistDeltaRDKR->Fill(dR0, ptD);
1038  }
1039  }
1040 }
1041 
1042 //_______________________________________________________________________________
1044 {
1045  // Fill some histogram with kinematic information of the D0 candidate.
1046 
1047  AliAODTrack *softPionTrack = static_cast<AliAODTrack*>(dstar->GetBachelor());
1048  Double_t ptD = dstar->Pt();
1049  Double_t aD = dstar->Phi();
1050  Double_t apis= softPionTrack->Phi();
1051 
1052  AliAODRecoDecayHF2Prong* D0fromDstar = dstar->Get2Prong();
1053  Double_t aD0 = D0fromDstar->Phi();
1054  Int_t isD0 = D0fromDstar->Charge()>0 ? kTRUE : kFALSE;
1055  Double_t aK = isD0 ? D0fromDstar->PhiProng(0) : D0fromDstar->PhiProng(1);
1056  Double_t api = isD0 ? D0fromDstar->PhiProng(1) : D0fromDstar->PhiProng(0);
1057  Double_t dRDD0 = DeltaR(dstar,D0fromDstar);
1058  Double_t dRDpis = DeltaR(dstar,softPionTrack);
1059  Double_t dRDpi = DeltaR(dstar, isD0 ? static_cast<AliVParticle*>(D0fromDstar->GetDaughter(1)) : static_cast<AliVParticle*>(D0fromDstar->GetDaughter(0)));
1060  Double_t dRDK = DeltaR(dstar, isD0 ? static_cast<AliVParticle*>(D0fromDstar->GetDaughter(0)) : static_cast<AliVParticle*>(D0fromDstar->GetDaughter(1)));
1061 
1062  if (isDstar) {
1063  fHistAlphaDDS ->Fill(aD-aD0, ptD);
1064  fHistAlphaDpisS->Fill(aD-apis, ptD);
1065  fHistAlphaDpiS ->Fill(aD-api, ptD);
1066  fHistAlphaDKS ->Fill(aD-aK, ptD);
1067 
1068  fHistDeltaRDDS ->Fill(dRDD0, ptD);
1069  fHistDeltaRDpisS->Fill(dRDpis, ptD);
1070  fHistDeltaRDpiS ->Fill(dRDpi, ptD);
1071  fHistDeltaRDKS ->Fill(dRDK, ptD);
1072  }
1073  else {
1074  fHistAlphaDDB ->Fill(aD-aD0, ptD);
1075  fHistAlphaDpisB->Fill(aD-apis, ptD);
1076  fHistAlphaDpiB ->Fill(aD-api, ptD);
1077  fHistAlphaDKB ->Fill(aD-aK, ptD);
1078 
1079  fHistDeltaRDDB ->Fill(dRDD0, ptD);
1080  fHistDeltaRDpisB->Fill(dRDpis, ptD);
1081  fHistDeltaRDpiB ->Fill(dRDpi, ptD);
1082  fHistDeltaRDKB ->Fill(dRDK, ptD);
1083  }
1084 }
1085 
1086 //_______________________________________________________________________________
1088 {
1089  // Fills the array of Dstar side band candidates.
1090 
1091  //select by track cuts the side band candidates (don't want mass cut)
1092  Int_t isSelected = fCuts->IsSelected(dstar, AliRDHFCuts::kTracks, fAodEvent);
1093  if (!isSelected) return;
1094 
1095  //add a reasonable cut on the invariant mass (e.g. (+-2\sigma, +-10 \sigma), with \sigma = fSigmaD0[bin])
1096  Int_t bin = fCuts->PtBin(dstar->Pt());
1097  if (bin < 0 || bin >= fCuts->GetNPtBins()) {
1098  AliError(Form("Pt %.3f out of bounds", dstar->Pt()));
1099  return;
1100  }
1101 
1102  const Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1103 
1104  //if data and Dstar from D0 side band
1105  if (((dstar->InvMassD0() < (mPDGD0-3.*fSigmaD0[bin])) && (dstar->InvMassD0() > (mPDGD0-10.*fSigmaD0[bin]))) /*left side band*/ ||
1106  ((dstar->InvMassD0() > (mPDGD0+3.*fSigmaD0[bin])) && (dstar->InvMassD0() < (mPDGD0+10.*fSigmaD0[bin]))) /*right side band*/) {
1107 
1108  if(!fMultCand) new ((*fSideBandArray)[fNSBCand]) AliAODRecoCascadeHF(*dstar);
1109  else if(fNSBCand==fAnalyseCand) new ((*fSideBandArray)[0]) AliAODRecoCascadeHF(*dstar);
1110  fNSBCand++;
1111  }
1112 }
1113 
1114 //_______________________________________________________________________________
1116 {
1117  // Set the mass limits using the PDG code and the given range.
1118 
1119  Float_t mass = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
1120 
1121  // compute the Delta mass for the D*
1122  if (fCandidateType==kDstartoKpipi) mass -= TDatabasePDG::Instance()->GetParticle(421)->Mass();
1123 
1124  fMinMass = mass - range;
1125  fMaxMass = mass + range;
1126 
1127  AliInfo(Form("Setting mass limits to %f, %f", fMinMass, fMaxMass));
1128  if ((fMinMass<0.) || (fMaxMass<=0.) || (fMaxMass<fMinMass)) {
1129  AliError(Form("Wrong mass limits! Task '%s' will be disabled!", GetName()));
1130  fInhibitTask = kTRUE;
1131  }
1132 }
1133 
1134 //_______________________________________________________________________________
1136 {
1137  // Set the mass limits.
1138 
1139  if (uplimit>lowlimit) {
1140  fMinMass = lowlimit;
1141  fMaxMass = uplimit;
1142  }
1143  else {
1144  printf("Error! Lower limit larger than upper limit!\n");
1145  fMinMass = uplimit - uplimit*0.2;
1146  fMaxMass = uplimit;
1147  }
1148 }
1149 
1150 //_______________________________________________________________________________
1152 {
1153  // Set the D0 width for the D* analysis.
1154 
1155  if (nptbins > 30) {
1156  AliWarning("Maximum number of bins allowed is 30!");
1157  return kFALSE;
1158  }
1159 
1160  if (!width) return kFALSE;
1161  for (Int_t ipt=0; ipt < nptbins; ipt++) fSigmaD0[ipt] = width[ipt];
1162 
1163  return kTRUE;
1164 }
1165 
1166 //_______________________________________________________________________________
1168 {
1169  // Allocate the histograms for the analysis.
1170 
1171  // Statistics
1172  fHistStat = new TH1I("fHistStat", "Statistics", 9, -0.5, 8.5);
1173  fHistStat->GetXaxis()->SetBinLabel(1, "N ev anal");
1174  fHistStat->GetXaxis()->SetBinLabel(2, "N ev sel");
1175  if(fUseMCInfo && fBuildRMEff){
1176  fHistStat->GetXaxis()->SetBinLabel(3, "N Gen D");
1177  fHistStat->GetXaxis()->SetBinLabel(4, "N Gen Sel D");
1178  fHistStat->GetXaxis()->SetBinLabel(5, "N D");
1179  fHistStat->GetXaxis()->SetBinLabel(6, "N cand sel cuts");
1180  fHistStat->GetXaxis()->SetBinLabel(7, "N rej no quark");
1181  fHistStat->GetXaxis()->SetBinLabel(8, "N rej from B");
1182  fHistStat->GetXaxis()->SetBinLabel(9, "N rej from D");
1183  }
1184  else {
1185  if (fUseMCInfo) {
1186  if (fUseReco) {
1187  fHistStat->GetXaxis()->SetBinLabel(3, "N D");
1188  }
1189  else {
1190  fHistStat->GetXaxis()->SetBinLabel(3, "N Gen D");
1191  }
1192  }
1193  else {
1194  fHistStat->GetXaxis()->SetBinLabel(3, "N Cand");
1195  }
1196  fHistStat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
1197  if (fCandidateType == kDstartoKpipi) {
1198  fHistStat->GetXaxis()->SetBinLabel(5, "N side band cand");
1199  }
1200  if (fUseMCInfo) {
1201  fHistStat->GetXaxis()->SetBinLabel(6, "N Background");
1202  fHistStat->GetXaxis()->SetBinLabel(7, "N rej no quark");
1203  fHistStat->GetXaxis()->SetBinLabel(8, "N rej from B");
1204  fHistStat->GetXaxis()->SetBinLabel(9, "N rej from D");
1205  }
1206  }
1207 
1208  fHistStat->SetNdivisions(1);
1209  fOutput->Add(fHistStat);
1210 
1211  fHistNCandEv = new TH1F("fHistNCandEv", "Number of candidates per event (after cuts);# cand/ev", 100, 0., 100.);
1212  fOutput->Add(fHistNCandEv);
1213 
1214  if ((fCandidateType == kDstartoKpipi) || fUseMCInfo) {
1215  fHistNSBCandEv = new TH1F("hnSBCandEv", "Number of side bands candidates per event (after cuts);# cand/ev", 100, 0.,100.);
1216  fOutput->Add(fHistNSBCandEv);
1217  }
1218 
1219  if (fCandidateType == kDstartoKpipi) {
1220  fHistPtPion = new TH1F("fHistPtPion", "Primary pions candidates pt", 500, 0., 10.);
1221  fHistPtPion->SetStats(kTRUE);
1222  fHistPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1223  fHistPtPion->GetYaxis()->SetTitle("entries");
1224  fOutput->Add(fHistPtPion);
1225  }
1226 
1227  // Invariant mass related histograms
1228  const Int_t nbinsmass = 200;
1229  const Int_t ptbinsD = 100;
1230  Float_t ptmin =0.;
1231  Float_t ptmax =50.;
1232  fHistInvMassPtD = new TH2F("fHistInvMassPtD", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, ptbinsD, ptmin, ptmax);
1233  fHistInvMassPtD->SetStats(kTRUE);
1234  fHistInvMassPtD->GetXaxis()->SetTitle("mass (GeV/c)");
1235  fHistInvMassPtD->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1236  fOutput->Add(fHistInvMassPtD);
1237 
1238  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
1239  fOutput->Add(fHistImpParS);
1240 
1241  fHistInvMassS = new TH1F("fHistInvMassS", "D invariant mass distribution (filled with fCandidateArray)", nbinsmass, fMinMass, fMaxMass);
1242  fHistInvMassS->SetStats(kTRUE);
1243  fHistInvMassS->GetXaxis()->SetTitle("mass (GeV/c)");
1244  fHistInvMassS->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1245  fOutput->Add(fHistInvMassS);
1246 
1247  const Int_t nbinsalpha = 200;
1248  Float_t minalpha = -TMath::Pi();
1249  Float_t maxalpha = TMath::Pi();
1250  const Int_t nbinsdeltaR = 200;
1251  Float_t mindeltaR = 0.;
1252  Float_t maxdeltaR = 10.;
1253 
1254 
1255 
1256  if (fUseMCInfo) {
1257  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
1258  fOutput->Add(fHistImpParB);
1259 
1260  fHistInvMassB = new TH1F("fHistInvMassB", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass);
1261  fHistInvMassB->SetStats(kTRUE);
1262  fHistInvMassB->GetXaxis()->SetTitle("mass (GeV/c)");
1263  fHistInvMassB->GetYaxis()->SetTitle("p_{T} (GeV/c)");
1264  fOutput->Add(fHistInvMassB);
1265 
1266  if (fCandidateType == kDstartoKpipi) {
1267 
1268  fHistAlphaDDS = new TH2F("fHistAlphaDDS","Angle D^{*}-D^{0} (Signal);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1269  fHistAlphaDpisS = new TH2F("fHistAlphaDpisS","Angle D^{*}-#pi_{soft} (Signal);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1270  fHistAlphaDpiS = new TH2F("fHistAlphaDpiS","Angle D^{*}-#pi (Signal);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1271  fHistAlphaDKS = new TH2F("fHistAlphaDKS","Angle D^{*}-K (Signal);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1272 
1273  fHistAlphaDDB = new TH2F("fHistAlphaDDB","Angle D^{*}-D^{0} (Background);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1274  fHistAlphaDpisB = new TH2F("fHistAlphaDpisB","Angle D^{*}-#pi_{soft} (Background);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1275  fHistAlphaDpiB = new TH2F("fHistAlphaDpiB","Angle D^{*}-#pi (Background);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1276  fHistAlphaDKB = new TH2F("fHistAlphaDKB","Angle D^{*}-K (Background);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1277 
1278  fHistDeltaRDDS = new TH2F("fHistDeltaRDDS","Angle D^{*}-D^{0} (Signal);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1279  fHistDeltaRDpisS = new TH2F("fHistDeltaRDpisS","Angle D^{*}-#pi_{soft} (Signal);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1280  fHistDeltaRDpiS = new TH2F("fHistDeltaRDpiS","Angle D^{*}-#pi (Signal);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1281  fHistDeltaRDKS = new TH2F("fHistDeltaRDKS","Angle D^{*}-K (Signal);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1282 
1283  fHistDeltaRDDB = new TH2F("fHistDeltaRDDB","Angle D^{*}-D^{0} (Background);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1284  fHistDeltaRDpisB = new TH2F("fHistDeltaRDpisB","Angle D^{*}-#pi_{soft} (Background);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1285  fHistDeltaRDpiB = new TH2F("fHistDeltaRDpiB","Angle D^{*}-#pi (Background);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1286  fHistDeltaRDKB = new TH2F("fHistDeltaRDKB","Angle D^{*}-K (Background);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1287 
1288  fOutput->Add(fHistAlphaDDS);
1289  fOutput->Add(fHistAlphaDpisS);
1290  fOutput->Add(fHistAlphaDpiS);
1291  fOutput->Add(fHistAlphaDKS);
1292  fOutput->Add(fHistAlphaDDB);
1293  fOutput->Add(fHistAlphaDpisB);
1294  fOutput->Add(fHistAlphaDpiB);
1295  fOutput->Add(fHistAlphaDKB);
1296 
1297  fOutput->Add(fHistDeltaRDDS);
1298  fOutput->Add(fHistDeltaRDpisS);
1299  fOutput->Add(fHistDeltaRDpiS);
1300  fOutput->Add(fHistDeltaRDKS);
1301  fOutput->Add(fHistDeltaRDDB);
1302  fOutput->Add(fHistDeltaRDpisB);
1303  fOutput->Add(fHistDeltaRDpiB);
1304  fOutput->Add(fHistDeltaRDKB);
1305  }
1306  else if (fCandidateType == kD0toKpi) {
1307 
1308  fHistAlphaDpiS = new TH2F("fHistAlphaDpiS","Angle D^{0}-#pi (Signal);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1309  fHistAlphaDKS = new TH2F("fHistAlphaDKS","Angle D^{0}-K (Signal);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1310  fHistAlphaDpiR = new TH2F("fHistAlphaDpiR","Angle D^{0}-#pi (Reflections);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1311  fHistAlphaDKR = new TH2F("fHistAlphaDKR","Angle D^{0}-K (Reflections);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1312 
1313  fHistAlphaDpiB = new TH2F("fHistAlphaDpiB","Angle D^{0}-#pi (Background);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1314  fHistAlphaDKB = new TH2F("fHistAlphaDKB","Angle D^{0}-K (Background);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
1315 
1316  fHistDeltaRDpiS = new TH2F("fHistDeltaRDpiS","Angle D^{0}-#pi (Signal);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1317  fHistDeltaRDKS = new TH2F("fHistDeltaRDKS","Angle D^{0}-K (Signal);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1318  fHistDeltaRDpiR = new TH2F("fHistDeltaRDpiR","Angle D^{0}-#pi (Reflections);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1319  fHistDeltaRDKR = new TH2F("fHistDeltaRDKR","Angle D^{0}-K (Reflections);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1320 
1321  fHistDeltaRDpiB = new TH2F("fHistDeltaRDpiB","Angle D^{0}-#pi (Background);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1322  fHistDeltaRDKB = new TH2F("fHistDeltaRDKB","Angle D^{0}-K (Background);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
1323 
1324  fOutput->Add(fHistAlphaDpiS);
1325  fOutput->Add(fHistAlphaDKS);
1326  fOutput->Add(fHistAlphaDpiR);
1327  fOutput->Add(fHistAlphaDKR);
1328  fOutput->Add(fHistAlphaDpiB);
1329  fOutput->Add(fHistAlphaDKB);
1330 
1331  fOutput->Add(fHistDeltaRDpiS);
1332  fOutput->Add(fHistDeltaRDKS);
1333  fOutput->Add(fHistDeltaRDpiR);
1334  fOutput->Add(fHistDeltaRDKR);
1335  fOutput->Add(fHistDeltaRDpiB);
1336  fOutput->Add(fHistDeltaRDKB);
1337  }
1338  }
1339 
1340  return kTRUE;
1341 }
1342 
1343 //_______________________________________________________________________________
1344 Float_t AliAnalysisTaskSEDmesonsFilterCJ::DeltaR(AliVParticle *p1, AliVParticle *p2) const
1345 {
1346  // Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1347 
1348  if (!p1 || !p2) return -1;
1349 
1350  Double_t phi1 = p1->Phi();
1351  Double_t eta1 = p1->Eta();
1352  Double_t phi2 = p2->Phi();
1353  Double_t eta2 = p2->Eta();
1354 
1355  Double_t dPhi = phi1 - phi2;
1356  if(dPhi < -(TMath::Pi())/2) dPhi = dPhi + TMath::TwoPi();
1357  if(dPhi > (3*(TMath::Pi()))/2) dPhi = dPhi - TMath::TwoPi();
1358 
1359  Double_t dEta = eta1 - eta2;
1360 
1361  Double_t deltaR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
1362 
1363  return deltaR;
1364 }
1365 
1366 //_______________________________________________________________________________
1368 {
1369  //
1370  // Add event tracks to a collection that already contains the D candidates, excluding the daughters of the D candidates
1371  //
1372 
1373  if (!tracks) return;
1374 
1375  TObjArray allDaughters(10);
1376  allDaughters.SetOwner(kFALSE);
1377 
1378  TIter next(coll);
1379  AliEmcalParticle* emcpart = 0;
1380  while ((emcpart = static_cast<AliEmcalParticle*>(next()))) {
1381  AliAODRecoDecay* reco = dynamic_cast<AliAODRecoDecay*>(emcpart->GetTrack());
1382  AliDebug(2, Form("Found a D meson candidtate with pT = %.3f, eta = %.3f, phi = %.3f", reco->Pt(), reco->Eta(), reco->Phi()));
1383  if (reco) AddDaughters(reco, allDaughters);
1384  }
1385 
1386  tracks->ResetCurrentID();
1387  AliVTrack* track = 0;
1388  Int_t n = coll->GetEntriesFast();
1389 
1390  while ((track = static_cast<AliVTrack*>(tracks->GetNextAcceptParticle()))) {
1391  if(fUseMCInfo && fUsePythia){
1392  AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(track);
1393  bool isInj = IsTrackInjected(aodtrack, fMCHeader, fMCarray);
1394  if(!isInj) continue;
1395  }
1396 
1397  if (allDaughters.Remove(track) == 0) {
1398  new ((*coll)[n]) AliEmcalParticle(track);
1399  n++;
1400  AliDebug(2, Form("Track %d (pT = %.3f, eta = %.3f, phi = %.3f) is included", tracks->GetCurrentID(), track->Pt(), track->Eta(), track->Phi()));
1401  }
1402  else {
1403  AliDebug(2, Form("Track %d (pT = %.3f, eta = %.3f, phi = %.3f) is excluded", tracks->GetCurrentID(), track->Pt(), track->Eta(), track->Phi()));
1404  }
1405  }
1406 }
1407 
1408 //_______________________________________________________________________________
1410 {
1411  // Add all the dauthers of cand in an array. Follows all the decay cascades.
1412 
1413  Int_t n = cand->GetNDaughters();
1414 
1415  //Printf("AddDaughters: the number of dauhters is %d", n);
1416 
1417  Int_t ntot = 0;
1418  Double_t pt = 0;
1419  for (Int_t i = 0; i < n; i++) {
1420  AliVTrack* track = dynamic_cast<AliVTrack*>(cand->GetDaughter(i));
1421  if (!track) continue;
1422 
1423  AliAODRecoDecay* cand2 = dynamic_cast<AliAODRecoDecay*>(track);
1424 
1425  if (cand2) {
1426  //Printf("Daughter pT = %.3f --> ", track->Pt());
1427  pt += AddDaughters(cand2, daughters);
1428  }
1429  else {
1430  if (!track->InheritsFrom("AliAODTrack")) {
1431  Printf("Warning: One of the daughters is not of type 'AliAODTrack' nor 'AliAODRecoDecay'.");
1432  continue;
1433  }
1434  //Printf("Daughter pT = %.3f", track->Pt());
1435  daughters.AddLast(track);
1436  pt += track->Pt();
1437  ntot++;
1438  }
1439  }
1440 
1441  //Printf("Total pt of the daughters = %.3f", pt);
1442 
1443  return pt;
1444 }
1445 //_______________________________________________________________________________
1447 {
1448  //
1449  // Add event tracks to a collection that already contains the D candidates, excluding the daughters of the D candidates
1450  //
1451 
1452  if (!mctracks) return;
1453 
1454  TObjArray allMCDaughters(10);
1455  allMCDaughters.SetOwner(kFALSE);
1456 
1457  AliAODMCParticle* mcD = (AliAODMCParticle*)coll->At(0);
1458 
1459  AddMCDaughters(mcD,allMCDaughters,fMCarray);
1460 
1461  mctracks->ResetCurrentID();
1462  AliAODMCParticle* mcpart = 0;
1463  Int_t n = coll->GetEntriesFast();
1464  while ((mcpart = static_cast<AliAODMCParticle*>(mctracks->GetNextAcceptParticle()))) {
1465  if(TMath::Abs(mcpart->Charge())==0) continue;
1466  if(fUsePythia){
1467  bool isInj = IsMCTrackInjected(mcpart, fMCHeader, fMCarray);
1468  if(!isInj) continue;
1469  }
1470  if (allMCDaughters.Remove(mcpart) == 0) {
1471  new ((*coll)[n]) AliAODMCParticle(*mcpart);
1472  n++;
1473  AliDebug(2, Form("Track %d (pT = %.3f, eta = %.3f, phi = %.3f) is included", mctracks->GetCurrentID(), mcpart->Pt(), mcpart->Eta(), mcpart->Phi()));
1474  }
1475  else {
1476  AliDebug(2, Form("Track %d (pT = %.3f, eta = %.3f, phi = %.3f) is excluded", mctracks->GetCurrentID(), mcpart->Pt(), mcpart->Eta(), mcpart->Phi()));
1477  }
1478  }
1479 }
1480 
1481 //_______________________________________________________________________________
1482 Double_t AliAnalysisTaskSEDmesonsFilterCJ::AddMCDaughters(AliAODMCParticle* mcDmeson, TObjArray& mcdaughters, TClonesArray* mcArray)
1483 {
1484  // Add all the dauthers of cand in an array. Follows all the decay cascades.
1485 
1486  Int_t n = mcDmeson->GetNDaughters();
1487 
1488  //Printf("AddDaughters: the number of dauhters is %d", n);
1489  Double_t pt = 0;
1490 
1491  for (Int_t i = 0; i < n; i++) {
1492  AliAODMCParticle* DDaughter = static_cast<AliAODMCParticle*>(mcArray->At(mcDmeson->GetDaughter(i)));
1493  if (!DDaughter) continue;
1494 
1495  if (DDaughter->GetNDaughters()>0) {
1496  //Printf("Daughter pT = %.3f --> ", track->Pt());
1497  pt += AddMCDaughters(DDaughter, mcdaughters, mcArray);
1498  }
1499  else {
1500  //Printf("Daughter pT = %.3f", track->Pt());
1501  mcdaughters.AddLast(DDaughter);
1502  pt += DDaughter->Pt();
1503  }
1504  }
1505 
1506  //Printf("Total pt of the daughters = %.3f", pt);
1507 
1508  return pt;
1509 }
1510 //_________________________________________________________________________________________________
1511 Int_t AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin(AliAODRecoDecay* cand, TClonesArray* mcArray)
1512 {
1513  // Checks whether the mother of the D meson candidate comes from a charm or a bottom quark.
1514 
1515  if (!mcArray) return -1;
1516 
1517  Int_t labDau0 = static_cast<AliVTrack*>(cand->GetDaughter(0))->GetLabel();
1518  if (labDau0 < 0) return -1;
1519 
1520  AliAODMCParticle* part = static_cast<AliAODMCParticle*>(mcArray->At(labDau0));
1521  return CheckOrigin(part, mcArray);
1522 }
1523 
1524 //_________________________________________________________________________________________________
1525 Int_t AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin(AliAODRecoDecay* cand, AliStack* stack)
1526 {
1527  // Checks whether the mother of the D meson candidate comes from a charm or a bottom quark.
1528 
1529  if (!stack) return -1;
1530 
1531  Int_t labDau0 = static_cast<AliVTrack*>(cand->GetDaughter(0))->GetLabel();
1532  if (labDau0 < 0) return -1;
1533 
1534  return CheckOrigin(labDau0, stack);
1535 }
1536 
1537 //_________________________________________________________________________________________________
1538 Int_t AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin(AliAODMCParticle* part, TClonesArray* mcArray)
1539 {
1540  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1541 
1542  if (!part) return -1;
1543  if (!mcArray) return -1;
1544 
1545  Int_t pdgGranma = 0;
1546  Int_t mother = part->GetMother();
1547  Int_t istep = 0;
1548  Int_t abspdgGranma = 0;
1549  Bool_t isFromB = kFALSE;
1550  Bool_t isQuarkFound = kFALSE;
1551 
1552  while (mother >= 0) {
1553  istep++;
1554  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1555  if (mcGranma >= 0) {
1556  pdgGranma = mcGranma->GetPdgCode();
1557  abspdgGranma = TMath::Abs(pdgGranma);
1558  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1559  isFromB = kTRUE;
1560  }
1561 
1562  if (abspdgGranma == 4 || abspdgGranma == 5) isQuarkFound = kTRUE;
1563  mother = mcGranma->GetMother();
1564  }
1565  else {
1566  ::Error("AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1567  break;
1568  }
1569  }
1570 
1571  if (isQuarkFound) {
1572  if (isFromB) {
1573  return kFromBottom;
1574  }
1575  else {
1576  return kFromCharm;
1577  }
1578  }
1579  else {
1580  return kQuarkNotFound;
1581  }
1582 }
1583 
1584 //_________________________________________________________________________________________________
1586 {
1587  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1588 
1589  if (!stack) return -1;
1590 
1591  TParticle* part = stack->Particle(ipart);
1592  if (!part) return -1;
1593 
1594  Int_t pdgGranma = 0;
1595  Int_t mother = part->GetFirstMother();
1596  Int_t istep = 0;
1597  Int_t abspdgGranma = 0;
1598  Bool_t isFromB = kFALSE;
1599  Bool_t isQuarkFound = kFALSE;
1600 
1601  while (mother >= 0) {
1602  istep++;
1603  TParticle* mcGranma = stack->Particle(mother);
1604  if (mcGranma >= 0) {
1605  pdgGranma = mcGranma->GetPdgCode();
1606  abspdgGranma = TMath::Abs(pdgGranma);
1607  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1608  isFromB = kTRUE;
1609  }
1610 
1611  if (abspdgGranma == 4 || abspdgGranma == 5) isQuarkFound = kTRUE;
1612  mother = mcGranma->GetFirstMother();
1613  }
1614  else {
1615  ::Error("AliAnalysisTaskSEDmesonsFilterCJ::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1616  break;
1617  }
1618  }
1619 
1620  if (isQuarkFound) {
1621  if (isFromB) {
1622  return kFromBottom;
1623  }
1624  else {
1625  return kFromCharm;
1626  }
1627  }
1628  else {
1629  return kQuarkNotFound;
1630  }
1631 }
1632 
1633 //_________________________________________________________________________________________________
1634 Int_t AliAnalysisTaskSEDmesonsFilterCJ::CheckDecayChannel(AliAODMCParticle* part, TClonesArray* mcArray)
1635 {
1636  // Determine the decay channel
1637 
1638  if (!part) return -1;
1639  if (!mcArray) return -1;
1640 
1642 
1643  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1644 
1645  if (part->GetNDaughters() == 2) {
1646 
1647  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1648  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1649 
1650  if (!d1 || !d2) {
1651  return decay;
1652  }
1653 
1654  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1655  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1656 
1657  if (absPdgPart == 421) { // D0 -> K pi
1658 
1659  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1660  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1661  decay = kDecayD0toKpi;
1662  }
1663  }
1664 
1665  if (absPdgPart == 413) { // D* -> D0 pi
1666 
1667  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1668  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1669  if (D0decay == kDecayD0toKpi) {
1670  decay = kDecayDStartoKpipi;
1671  }
1672  }
1673 
1674  if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1675  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1676  if (D0decay == kDecayD0toKpi) {
1677  decay = kDecayDStartoKpipi;
1678  }
1679  }
1680  }
1681  }
1682 
1683  return decay;
1684 }
1685 
1686 //_________________________________________________________________________________________________
1688 {
1689  // Determine the decay channel
1690 
1691  if (!stack) return -1;
1692 
1693  TParticle* part = stack->Particle(ipart);
1694 
1695  if (!part) return -1;
1696 
1698 
1699  if (part->GetNDaughters() == 2) {
1700 
1701  Int_t id1 = part->GetDaughter(0);
1702  Int_t id2 = part->GetDaughter(1);
1703 
1704  TParticle* d1 = stack->Particle(id1);
1705  TParticle* d2 = stack->Particle(id2);
1706 
1707  if (!d1 || !d2) {
1708  return decay;
1709  }
1710 
1711  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1712  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1713 
1714 
1715  if (part->GetPdgCode() == 421) { // D0 -> K pi
1716 
1717  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1718  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1719  decay = kDecayD0toKpi;
1720  }
1721  }
1722 
1723  if (part->GetPdgCode() == 413) { // D* -> D0 pi
1724 
1725  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1726  Int_t D0decay = CheckDecayChannel(id1, stack);
1727  if (D0decay == kDecayD0toKpi) {
1728  decay = kDecayDStartoKpipi;
1729  }
1730  }
1731 
1732  if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1733  Int_t D0decay = CheckDecayChannel(id2, stack);
1734  if (D0decay == kDecayD0toKpi) {
1735  decay = kDecayDStartoKpipi;
1736  }
1737  }
1738  }
1739  }
1740 
1741  return decay;
1742 
1743 }
1744 
1745 
1746 //_____________________________________________________________________
1747 void AliAnalysisTaskSEDmesonsFilterCJ::GetTrackPrimaryGenerator(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){
1748 
1750 
1751  Int_t lab=TMath::Abs(track->GetLabel());
1752  nameGen=AliVertexingHFUtils::GetGenerator(lab,header);
1753 
1754  // Int_t countControl=0;
1755 
1756  while(nameGen.IsWhitespace()){
1757  AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);
1758  if(!mcpart){
1759  printf("AliAnalysisTaskMultCheck::IsTrackInjected - BREAK: No valid AliAODMCParticle at label %i\n",lab);
1760  break;
1761  }
1762  Int_t mother = mcpart->GetMother();
1763  if(mother<0){
1764  printf("AliAnalysisTaskMultCheck::IsTrackInjected - BREAK: Reached primary particle without valid mother\n");
1765  break;
1766  }
1767  lab=mother;
1768  nameGen=AliVertexingHFUtils::GetGenerator(mother,header);
1769  // countControl++;
1770  // if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops
1771  // printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Protection from infinite loop active\n");
1772  // break;
1773  // }
1774  }
1775 
1776  return;
1777 }
1778 //----------------------------------------------------------------------
1779 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::IsTrackInjected(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC){
1781  TString nameGen;
1782  GetTrackPrimaryGenerator(track,header,arrayMC,nameGen);
1783 
1784  if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
1785 
1786  return kTRUE;
1787 }
1788 
1789 
1790 //_____________________________________________________________________
1791 void AliAnalysisTaskSEDmesonsFilterCJ::GetMCTrackPrimaryGenerator(AliAODMCParticle *track,AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){
1792 
1794 
1795  Int_t lab=TMath::Abs(track->GetLabel());
1796  nameGen=AliVertexingHFUtils::GetGenerator(lab,header);
1797 
1798  // Int_t countControl=0;
1799 
1800  while(nameGen.IsWhitespace()){
1801  AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);
1802  if(!mcpart){
1803  printf("AliAnalysisTaskMultCheck::IsTrackInjected - BREAK: No valid AliAODMCParticle at label %i\n",lab);
1804  break;
1805  }
1806  Int_t mother = mcpart->GetMother();
1807  if(mother<0){
1808  printf("AliAnalysisTaskMultCheck::IsTrackInjected - BREAK: Reached primary particle without valid mother\n");
1809  break;
1810  }
1811  lab=mother;
1812  nameGen=AliVertexingHFUtils::GetGenerator(mother,header);
1813  // countControl++;
1814  // if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops
1815  // printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Protection from infinite loop active\n");
1816  // break;
1817  // }
1818  }
1819 
1820  return;
1821 }
1822 //----------------------------------------------------------------------
1823 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::IsMCTrackInjected(AliAODMCParticle *track,AliAODMCHeader *header,TClonesArray *arrayMC){
1825  TString nameGen;
1826  GetMCTrackPrimaryGenerator(track,header,arrayMC,nameGen);
1827 
1828  if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
1829 
1830  return kTRUE;
1831 }
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
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.
ClassImp(AliAnalysisTaskSEDmesonsFilterCJ) AliAnalysisTaskSEDmesonsFilterCJ
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...