]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliAnalysisTaskSEDStarSpectra.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Add selection on track filter mask in D meson tasks
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDStarSpectra.cxx
CommitLineData
63147a75 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 * appeuear 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/* $Id$ */
17
18//
19//
20// Base class for DStar Analysis
21//
22//
23// The D* spectra study is done in pt bins:
24// [0,0.5] [0.5,1] [1,2] [2,3] [3,4] [4,5] [5,6] [6,7] [7,8],
25// [8,10],[10,12], [12,16], [16,20] and [20,24]
26//
27// Cuts arew centralized in AliRDHFCutsDStartoKpipi
28// Side Band and like sign background are implemented in the macro
29//
30//-----------------------------------------------------------------------
31//
32// Author A.Grelli
33// ERC-QGP Utrecht University - a.grelli@uu.nl,
34// Author Y.Wang
35// University of Heidelberg - yifei@physi.uni-heidelberg.de
36// Author C.Ivan
37// ERC-QGP Utrecht University - c.ivan@uu.nl,
38//
39//-----------------------------------------------------------------------
40
41#include <TSystem.h>
42#include <TParticle.h>
43#include <TH1I.h>
44#include "TROOT.h"
45#include <TDatabasePDG.h>
46#include <AliAnalysisDataSlot.h>
47#include <AliAnalysisDataContainer.h>
48#include "AliRDHFCutsDStartoKpipi.h"
49#include "AliStack.h"
50#include "AliMCEvent.h"
51#include "AliAnalysisManager.h"
52#include "AliAODMCHeader.h"
53#include "AliAODHandler.h"
54#include "AliLog.h"
55#include "AliAODVertex.h"
56#include "AliAODRecoDecay.h"
57#include "AliAODRecoDecayHF.h"
58#include "AliAODRecoCascadeHF.h"
59#include "AliAODRecoDecayHF2Prong.h"
60#include "AliAnalysisVertexingHF.h"
61#include "AliESDtrack.h"
62#include "AliAODMCParticle.h"
63#include "AliAnalysisTaskSE.h"
64#include "AliAnalysisTaskSEDStarSpectra.h"
65#include "AliNormalizationCounter.h"
66
67ClassImp(AliAnalysisTaskSEDStarSpectra)
68
69//__________________________________________________________________________
70AliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra():
71 AliAnalysisTaskSE(),
72 fEvents(0),
73 fAnalysis(0),
74 fD0Window(0),
75 fPeakWindow(0),
76 fUseMCInfo(kFALSE),
77 fDoSearch(kFALSE),
78 fOutput(0),
79 fOutputAll(0),
80 fOutputPID(0),
81 fNSigma(3),
82 fCuts(0),
83 fCEvents(0),
84 fTrueDiff2(0),
85 fDeltaMassD1(0),
86 fCounter(0),
87 fDoImpParDstar(kFALSE),
88 fNImpParBins(400),
89 fLowerImpPar(-2000.),
90 fHigherImpPar(2000.)
91{
92 //
93 // Default ctor
94 //
95}
96//___________________________________________________________________________
97AliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra(const Char_t* name, AliRDHFCutsDStartoKpipi* cuts) :
98 AliAnalysisTaskSE(name),
99 fEvents(0),
100 fAnalysis(0),
101 fD0Window(0),
102 fPeakWindow(0),
103 fUseMCInfo(kFALSE),
104 fDoSearch(kFALSE),
105 fOutput(0),
106 fOutputAll(0),
107 fOutputPID(0),
108 fNSigma(3),
109 fCuts(0),
110 fCEvents(0),
111 fTrueDiff2(0),
112 fDeltaMassD1(0),
113 fCounter(0),
114 fDoImpParDstar(kFALSE),
115 fNImpParBins(400),
116 fLowerImpPar(-2000.),
117 fHigherImpPar(2000.)
118{
119 //
120 // Constructor. Initialization of Inputs and Outputs
121 //
122 Info("AliAnalysisTaskSEDStarSpectra","Calling Constructor");
123
124 fCuts=cuts;
125
126 DefineOutput(1,TList::Class()); //conters
127 DefineOutput(2,TList::Class()); //All Entries output
128 DefineOutput(3,TList::Class()); //3sigma PID output
129 DefineOutput(4,AliRDHFCutsDStartoKpipi::Class()); //My private output
130 DefineOutput(5,AliNormalizationCounter::Class()); // normalization
131}
132
133//___________________________________________________________________________
134AliAnalysisTaskSEDStarSpectra::~AliAnalysisTaskSEDStarSpectra() {
135 //
136 // destructor
137 //
138 Info("~AliAnalysisTaskSEDStarSpectra","Calling Destructor");
139
140 if (fOutput) {
141 delete fOutput;
142 fOutput = 0;
143 }
144 if (fOutputAll) {
145 delete fOutputAll;
146 fOutputAll = 0;
147 }
148 if (fOutputPID) {
149 delete fOutputPID;
150 fOutputPID = 0;
151 }
152 if (fCuts) {
153 delete fCuts;
154 fCuts = 0;
155 }
156 if(fCEvents){
157 delete fCEvents;
158 fCEvents =0;
159 }
160 if(fDeltaMassD1){
161 delete fDeltaMassD1;
162 fDeltaMassD1 =0;
163 }
164 for(Int_t i=0; i<5; i++){
165 delete fHistMassPtImpParTCDs[i];
166 }
167}
168//_________________________________________________
169void AliAnalysisTaskSEDStarSpectra::Init(){
170 //
171 // Initialization
172 //
173
174 if(fDebug > 1) printf("AnalysisTaskSEDStarSpectra::Init() \n");
175 AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
176 // Post the data
177 PostData(4,copyfCuts);
178
179 return;
180}
181
182//_________________________________________________
183void AliAnalysisTaskSEDStarSpectra::UserExec(Option_t *)
184{
185 // user exec
186 if (!fInputEvent) {
187 Error("UserExec","NO EVENT FOUND!");
188 return;
189 }
190
191 fEvents++;
192
193 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
194 TClonesArray *arrayDStartoD0pi=0;
195
196 fCEvents->Fill(1);
197
198 if(!aodEvent && AODEvent() && IsStandardAOD()) {
199 // In case there is an AOD handler writing a standard AOD, use the AOD
200 // event in memory rather than the input (ESD) event.
201 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
202 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
203 // have to taken from the AOD event hold by the AliAODExtension
204 AliAODHandler* aodHandler = (AliAODHandler*)
205 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
206 if(aodHandler->GetExtensions()) {
207 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
208 AliAODEvent *aodFromExt = ext->GetAOD();
209 arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
210 }
211 } else {
212 arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
213 }
214
215 // fix for temporary bug in ESDfilter
216 // the AODs with null vertex pointer didn't pass the PhysSel
217 if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
218 fCEvents->Fill(2);
219
220 fCounter->StoreEvent(aodEvent,fCuts,fUseMCInfo);
221
222 // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
223 TString trigclass=aodEvent->GetFiredTriggerClasses();
224 if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")||trigclass.Contains("C0SMH-B-NOPF-ALL")) fCEvents->Fill(5);
225
226 if(!fCuts->IsEventSelected(aodEvent)) {
227 if(fCuts->GetWhyRejection()==6) // rejected for Z vertex
228 fCEvents->Fill(6);
229 return;
230 }
231
232 Bool_t isEvSel=fCuts->IsEventSelected(aodEvent);
233 if(!isEvSel) return;
234
235 // Load the event
236 // AliInfo(Form("Event %d",fEvents));
237 //if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
238
239 // counters for efficiencies
240 Int_t icountReco = 0;
241
242 //D* and D0 prongs needed to MatchToMC method
243 Int_t pdgDgDStartoD0pi[2]={421,211};
244 Int_t pdgDgD0toKpi[2]={321,211};
245
246 // AOD primary vertex
247 AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
248 if(!vtx1) return;
249 if(vtx1->GetNContributors()<1) return;
250
251 if (!arrayDStartoD0pi){
252 AliInfo("Could not find array of HF vertices, skipping the event");
253 return;
254 }else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
255
256 Int_t nSelectedAna =0;
257 Int_t nSelectedProd =0;
258
259 // loop over the tracks to search for candidates soft pion
260 for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
261
262 // D* candidates and D0 from D*
263 AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
264 if(!dstarD0pi->GetSecondaryVtx()) continue;
265 AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
266 if (!theD0particle) continue;
267
268 Int_t isDStar = 0;
269 TClonesArray *mcArray = 0;
270 AliAODMCHeader *mcHeader=0;
271
272 Bool_t isPrimary=kTRUE;
273 Float_t truePt=0.;
274 Float_t xDecay=0.;
275 Float_t yDecay=0.;
276 Float_t zDecay=0.;
277 Float_t pdgCode=-2;
278 Float_t trueImpParXY=0.;
279
280 // mc analysis
281 if(fUseMCInfo){
282 //MC array need for maching
283 mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
284 if (!mcArray) {
285 AliError("Could not find Monte-Carlo in AOD");
286 return;
287 }
288 // load MC header
289 mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
290 if(!mcHeader) {
291 printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
292 return;
293 }
294 // find associated MC particle for D* ->D0toKpi
295 Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
296 if(mcLabel>=0){
297
298 AliAODMCParticle *partDSt = (AliAODMCParticle*)mcArray->At(mcLabel);
299 Int_t checkOrigin = CheckOrigin(mcArray,partDSt);
300 if(checkOrigin==5) isPrimary=kFALSE;
301 AliAODMCParticle *dg0 = (AliAODMCParticle*)mcArray->At(partDSt->GetDaughter(0));
302 AliAODMCParticle *dg0_1 = (AliAODMCParticle*)mcArray->At(dg0->GetDaughter(0));
303 truePt=dg0->Pt();
304 xDecay=dg0_1->Xv();
305 yDecay=dg0_1->Yv();
306 zDecay=dg0_1->Zv();
307 pdgCode=TMath::Abs(partDSt->GetPdgCode());
308 if(!isPrimary){
309 trueImpParXY=GetTrueImpactParameterD0(mcHeader,mcArray,dg0)*1000.;
310 }
311 isDStar = 1;
312 }else{
313 pdgCode=-1;
314 }
315 }
316
317 Int_t ptbin=fCuts->PtBin(dstarD0pi->Pt());
318
319 // quality selction on tracks and region of interest
320 Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
321 if(!isTkSelected) continue;
322
323 if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
324
325
326 //histos for impact par studies - D0!!!
327 Double_t ptCand = dstarD0pi->Get2Prong()->Pt();
328 Double_t invMass=dstarD0pi->InvMassD0();
329 Double_t impparXY=dstarD0pi->Get2Prong()->ImpParXY()*10000.;
330
331 Double_t arrayForSparse[3]={invMass,ptCand,impparXY};
332 Double_t arrayForSparseTrue[3]={invMass,ptCand,trueImpParXY};
333
334 // set the D0 search window bin by bin - useful to calculate side band bkg
335 if (ptbin==0){
336 if(fAnalysis==1){
337 fD0Window=0.035;
338 fPeakWindow=0.03;
339 }else{
340 fD0Window=0.020;
341 fPeakWindow=0.0018;
342 }
343 }
344 if (ptbin==1){
345 if(fAnalysis==1){
346 fD0Window=0.035;
347 fPeakWindow=0.03;
348 }else{
349 fD0Window=0.020;
350 fPeakWindow=0.0018;
351 }
352 }
353 if (ptbin==2){
354 if(fAnalysis==1){
355 fD0Window=0.035;
356 fPeakWindow=0.03;
357 }else{
358 fD0Window=0.020;
359 fPeakWindow=0.0018;
360 }
361 }
362 if (ptbin==3){
363 if(fAnalysis==1){
364 fD0Window=0.035;
365 fPeakWindow=0.03;
366 }else{
367 fD0Window=0.022;
368 fPeakWindow=0.0016;
369 }
370 }
371 if (ptbin==4){
372 if(fAnalysis==1){
373 fD0Window=0.035;
374 fPeakWindow=0.03;
375 }else{
376 fD0Window=0.026;
377 fPeakWindow=0.0014;
378 }
379 }
380 if (ptbin==5){
381 if(fAnalysis==1){
382 fD0Window=0.045;
383 fPeakWindow=0.03;
384 }else{
385 fD0Window=0.026;
386 fPeakWindow=0.0014;
387 }
388 }
389 if (ptbin==6){
390 if(fAnalysis==1){
391 fD0Window=0.045;
392 fPeakWindow=0.03;
393 }else{
394 fD0Window=0.026;
395 fPeakWindow=0.006;
396 }
397 }
398 if (ptbin==7){
399 if(fAnalysis==1){
400 fD0Window=0.055;
401 fPeakWindow=0.03;
402 }else{
403 fD0Window=0.026;
404 fPeakWindow=0.006;
405 }
406 }
407 if (ptbin>7){
408 if(fAnalysis==1){
409 fD0Window=0.074;
410 fPeakWindow=0.03;
411 }else{
412 fD0Window=0.026;
413 fPeakWindow=0.006;
414 }
415 }
416
417 nSelectedProd++;
418 nSelectedAna++;
419
420 // check that we are close to signal in the DeltaM - here to save time for PbPb
421 Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
422 Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
423 Double_t invmassDelta = dstarD0pi->DeltaInvMass();
424
425 if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))>fPeakWindow) continue;
426
427
428 Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
429
430 // after cuts
431 if(fDoImpParDstar && isSelected){
432 fHistMassPtImpParTCDs[0]->Fill(arrayForSparse);
433 if(isPrimary) fHistMassPtImpParTCDs[1]->Fill(arrayForSparse);
434 else{
435 fHistMassPtImpParTCDs[2]->Fill(arrayForSparse);
436 fHistMassPtImpParTCDs[3]->Fill(arrayForSparseTrue);
437 }
438 }
439
440 // fill PID
441 FillSpectrum(dstarD0pi,isDStar,fCuts,isSelected,fOutputPID);
442 SideBandBackground(dstarD0pi,fCuts,isSelected, fOutputPID);
443 //WrongSignForDStar(dstarD0pi,fCuts,fOutputPID);
444
445 //swich off the PID selection
446 fCuts->SetUsePID(kFALSE);
447 Int_t isSelectedNoPID=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
448 fCuts->SetUsePID(kTRUE);
449
450 FillSpectrum(dstarD0pi,isDStar,fCuts,isSelectedNoPID,fOutputPID);
451 SideBandBackground(dstarD0pi,fCuts,isSelectedNoPID, fOutputPID);
452
453 // rare D search ------
454 if(fDoSearch){
455 TLorentzVector LorentzTrack1(0,0,0,0); // lorentz 4 vector
456 TLorentzVector LorentzTrack2(0,0,0,0); // lorentz 4 vector
457
458 for (Int_t i=0; i<aodEvent->GetNTracks(); i++){
459
460 AliAODTrack* aodTrack = aodEvent->GetTrack(i);
461
462 if(dstarD0pi->Charge() == aodTrack->Charge()) continue;
463 if((!(aodTrack->GetStatus()&AliESDtrack::kITSrefit)|| (!(aodTrack->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
464 if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))>0.02) continue;
465
466 //build the D1 mass
467 Double_t mass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
468
469 LorentzTrack1.SetPxPyPzE( dstarD0pi->Px(),dstarD0pi->Py(), dstarD0pi->Pz(), dstarD0pi->E(413) );
470 LorentzTrack2.SetPxPyPzE( aodTrack->Px(),aodTrack->Py(), aodTrack->Pz(),aodTrack->E(mass) );
471
472 //D1 mass
473 Double_t d1mass = ((LorentzTrack1+LorentzTrack2).M());
474 //mass difference - at 0.4117 and 0.4566
475 fDeltaMassD1->Fill(d1mass-dstarD0pi->InvMassDstarKpipi());
476 }
477 }
478
479 if(isDStar == 1) {
480 fTrueDiff2->Fill(dstarD0pi->Pt(),dstarD0pi->DeltaInvMass());
481 }
482
483 }
484
485 fCounter->StoreCandidates(aodEvent,nSelectedProd,kTRUE);
486 fCounter->StoreCandidates(aodEvent,nSelectedAna,kFALSE);
487
488 AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
489
490 PostData(1,fOutput);
491 PostData(2,fOutputAll);
492 PostData(3,fOutputPID);
493 PostData(5,fCounter);
494
495}
496//________________________________________ terminate ___________________________
497void AliAnalysisTaskSEDStarSpectra::Terminate(Option_t*)
498{
499 // The Terminate() function is the last function to be called during
500 // a query. It always runs on the client, it can be used to present
501 // the results graphically or save the results to file.
502
503 //Info("Terminate","");
504 AliAnalysisTaskSE::Terminate();
505
506 fOutput = dynamic_cast<TList*> (GetOutputData(1));
507 if (!fOutput) {
508 printf("ERROR: fOutput not available\n");
509 return;
510 }
511
512 fCEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fCEvents"));
513 fDeltaMassD1 = dynamic_cast<TH1F*>(fOutput->FindObject("fDeltaMassD1"));
514 fTrueDiff2 = dynamic_cast<TH2F*>(fOutput->FindObject("fTrueDiff2"));
515
516 fOutputAll = dynamic_cast<TList*> (GetOutputData(1));
517 if (!fOutputAll) {
518 printf("ERROR: fOutputAll not available\n");
519 return;
520 }
521 fOutputPID = dynamic_cast<TList*> (GetOutputData(2));
522 if (!fOutputPID) {
523 printf("ERROR: fOutputPID not available\n");
524 return;
525 }
526
527
528 return;
529}
530//___________________________________________________________________________
531void AliAnalysisTaskSEDStarSpectra::UserCreateOutputObjects() {
532 // output
533 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
534
535 //slot #1
536 //OpenFile(1);
537 fOutput = new TList();
538 fOutput->SetOwner();
539 fOutput->SetName("chist0");
540
541 fOutputAll = new TList();
542 fOutputAll->SetOwner();
543 fOutputAll->SetName("listAll");
544
545 fOutputPID = new TList();
546 fOutputPID->SetOwner();
547 fOutputPID->SetName("listPID");
548
549 // define histograms
550 DefineHistograms();
551
552 //Counter for Normalization
553 fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));
554
555 if(fDoImpParDstar) CreateImpactParameterHistos();
556
557 PostData(1,fOutput);
558 PostData(2,fOutputAll);
559 PostData(3,fOutputPID);
560
561 return;
562}
563//___________________________________ hiostograms _______________________________________
564void AliAnalysisTaskSEDStarSpectra::DefineHistograms(){
565
566 fCEvents = new TH1F("fCEvents","conter",11,0,11);
567 fCEvents->SetStats(kTRUE);
568 fCEvents->GetXaxis()->SetTitle("1");
569 fCEvents->GetYaxis()->SetTitle("counts");
570 fOutput->Add(fCEvents);
571
572 fTrueDiff2 = new TH2F("DiffDstar_pt","True Reco diff vs pt",200,0,15,900,0,0.3);
573 fOutput->Add(fTrueDiff2);
574
575 fDeltaMassD1 = new TH1F("DeltaMassD1","delta mass d1",600,0,0.8);
576 fOutput->Add(fDeltaMassD1);
577
578 const Int_t nhist=14;
579 TString nameMass=" ", nameSgn=" ", nameBkg=" ";
580
581 for(Int_t i=-2;i<nhist;i++){
582 nameMass="histDeltaMass_";
583 nameMass+=i+1;
584 nameSgn="histDeltaSgn_";
585 nameSgn+=i+1;
586 nameBkg="histDeltaBkg_";
587 nameBkg+=i+1;
588
589 if (i==-2) {
590 nameMass="histDeltaMass";
591 nameSgn="histDeltaSgn";
592 nameBkg="histDeltaBkg";
593 }
594
595 TH1F* spectrumMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} invariant mass; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
596 TH1F* spectrumSgn = new TH1F(nameSgn.Data(), "D^{*}-D^{0} Signal invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
597 TH1F* spectrumBkg = new TH1F(nameBkg.Data(), "D^{*}-D^{0} Background invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
598
599 nameMass="histD0Mass_";
600 nameMass+=i+1;
601 nameSgn="histD0Sgn_";
602 nameSgn+=i+1;
603 nameBkg="histD0Bkg_";
604 nameBkg+=i+1;
605
606 if (i==-2) {
607 nameMass="histD0Mass";
608 nameSgn="histD0Sgn";
609 nameBkg="histD0Bkg";
610 }
611
612 TH1F* spectrumD0Mass = new TH1F(nameMass.Data(),"D^{0} invariant mass; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
613 TH1F* spectrumD0Sgn = new TH1F(nameSgn.Data(), "D^{0} Signal invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
614 TH1F* spectrumD0Bkg = new TH1F(nameBkg.Data(), "D^{0} Background invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
615
616 nameMass="histDstarMass_";
617 nameMass+=i+1;
618 nameSgn="histDstarSgn_";
619 nameSgn+=i+1;
620 nameBkg="histDstarBkg_";
621 nameBkg+=i+1;
622
623 if (i==-2) {
624 nameMass="histDstarMass";
625 nameSgn="histDstarSgn";
626 nameBkg="histDstarBkg";
627 }
628
629 TH1F* spectrumDstarMass = new TH1F(nameMass.Data(),"D^{*} invariant mass; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
630 TH1F* spectrumDstarSgn = new TH1F(nameSgn.Data(), "D^{*} Signal invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
631 TH1F* spectrumDstarBkg = new TH1F(nameBkg.Data(), "D^{*} Background invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
632
633 nameMass="histSideBandMass_";
634 nameMass+=i+1;
635 if (i==-2) {
636 nameMass="histSideBandMass";
637 }
638
639 TH1F* spectrumSideBandMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} sideband mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
640
641 nameMass="histWrongSignMass_";
642 nameMass+=i+1;
643 if (i==-2) {
644 nameMass="histWrongSignMass";
645 }
646
647 TH1F* spectrumWrongSignMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} wrongsign mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
648
649
650 spectrumMass->Sumw2();
651 spectrumSgn->Sumw2();
652 spectrumBkg->Sumw2();
653
654 spectrumMass->SetLineColor(6);
655 spectrumSgn->SetLineColor(2);
656 spectrumBkg->SetLineColor(4);
657
658 spectrumMass->SetMarkerStyle(20);
659 spectrumSgn->SetMarkerStyle(20);
660 spectrumBkg->SetMarkerStyle(20);
661 spectrumMass->SetMarkerSize(0.6);
662 spectrumSgn->SetMarkerSize(0.6);
663 spectrumBkg->SetMarkerSize(0.6);
664 spectrumMass->SetMarkerColor(6);
665 spectrumSgn->SetMarkerColor(2);
666 spectrumBkg->SetMarkerColor(4);
667
668 spectrumD0Mass->Sumw2();
669 spectrumD0Sgn->Sumw2();
670 spectrumD0Bkg->Sumw2();
671
672 spectrumD0Mass->SetLineColor(6);
673 spectrumD0Sgn->SetLineColor(2);
674 spectrumD0Bkg->SetLineColor(4);
675
676 spectrumD0Mass->SetMarkerStyle(20);
677 spectrumD0Sgn->SetMarkerStyle(20);
678 spectrumD0Bkg->SetMarkerStyle(20);
679 spectrumD0Mass->SetMarkerSize(0.6);
680 spectrumD0Sgn->SetMarkerSize(0.6);
681 spectrumD0Bkg->SetMarkerSize(0.6);
682 spectrumD0Mass->SetMarkerColor(6);
683 spectrumD0Sgn->SetMarkerColor(2);
684 spectrumD0Bkg->SetMarkerColor(4);
685
686 spectrumDstarMass->Sumw2();
687 spectrumDstarSgn->Sumw2();
688 spectrumDstarBkg->Sumw2();
689
690 spectrumDstarMass->SetLineColor(6);
691 spectrumDstarSgn->SetLineColor(2);
692 spectrumDstarBkg->SetLineColor(4);
693
694 spectrumDstarMass->SetMarkerStyle(20);
695 spectrumDstarSgn->SetMarkerStyle(20);
696 spectrumDstarBkg->SetMarkerStyle(20);
697 spectrumDstarMass->SetMarkerSize(0.6);
698 spectrumDstarSgn->SetMarkerSize(0.6);
699 spectrumDstarBkg->SetMarkerSize(0.6);
700 spectrumDstarMass->SetMarkerColor(6);
701 spectrumDstarSgn->SetMarkerColor(2);
702 spectrumDstarBkg->SetMarkerColor(4);
703
704 spectrumSideBandMass->Sumw2();
705 spectrumSideBandMass->SetLineColor(4);
706 spectrumSideBandMass->SetMarkerStyle(20);
707 spectrumSideBandMass->SetMarkerSize(0.6);
708 spectrumSideBandMass->SetMarkerColor(4);
709
710 spectrumWrongSignMass->Sumw2();
711 spectrumWrongSignMass->SetLineColor(4);
712 spectrumWrongSignMass->SetMarkerStyle(20);
713 spectrumWrongSignMass->SetMarkerSize(0.6);
714 spectrumWrongSignMass->SetMarkerColor(4);
715
716 TH1F* allMass = (TH1F*)spectrumMass->Clone();
717 TH1F* allSgn = (TH1F*)spectrumSgn->Clone();
718 TH1F* allBkg = (TH1F*)spectrumBkg->Clone();
719
720 TH1F* pidMass = (TH1F*)spectrumMass->Clone();
721 TH1F* pidSgn = (TH1F*)spectrumSgn->Clone();
722 TH1F* pidBkg = (TH1F*)spectrumBkg->Clone();
723
724 fOutputAll->Add(allMass);
725 fOutputAll->Add(allSgn);
726 fOutputAll->Add(allBkg);
727
728 fOutputPID->Add(pidMass);
729 fOutputPID->Add(pidSgn);
730 fOutputPID->Add(pidBkg);
731
732 TH1F* allD0Mass = (TH1F*)spectrumD0Mass->Clone();
733 TH1F* allD0Sgn = (TH1F*)spectrumD0Sgn->Clone();
734 TH1F* allD0Bkg = (TH1F*)spectrumD0Bkg->Clone();
735
736 TH1F* pidD0Mass = (TH1F*)spectrumD0Mass->Clone();
737 TH1F* pidD0Sgn = (TH1F*)spectrumD0Sgn->Clone();
738 TH1F* pidD0Bkg = (TH1F*)spectrumD0Bkg->Clone();
739
740 fOutputAll->Add(allD0Mass);
741 fOutputAll->Add(allD0Sgn);
742 fOutputAll->Add(allD0Bkg);
743
744 fOutputPID->Add(pidD0Mass);
745 fOutputPID->Add(pidD0Sgn);
746 fOutputPID->Add(pidD0Bkg);
747
748 TH1F* allDstarMass = (TH1F*)spectrumDstarMass->Clone();
749 TH1F* allDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
750 TH1F* allDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
751
752 TH1F* pidDstarMass = (TH1F*)spectrumDstarMass->Clone();
753 TH1F* pidDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
754 TH1F* pidDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
755
756 fOutputAll->Add(allDstarMass);
757 fOutputAll->Add(allDstarSgn);
758 fOutputAll->Add(allDstarBkg);
759
760 fOutputPID->Add(pidDstarMass);
761 fOutputPID->Add(pidDstarSgn);
762 fOutputPID->Add(pidDstarBkg);
763
764 TH1F* allSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
765 TH1F* pidSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
766
767 fOutputAll->Add(allSideBandMass);
768 fOutputPID->Add(pidSideBandMass);
769
770 TH1F* allWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
771 TH1F* pidWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
772
773 fOutputAll->Add(allWrongSignMass);
774 fOutputPID->Add(pidWrongSignMass);
775
776 }
777
778 // pt spectra
779 nameMass="ptMass";
780 nameSgn="ptSgn";
781 nameBkg="ptBkg";
782
783 TH1F* ptspectrumMass = new TH1F(nameMass.Data(),"D^{*} p_{T}; p_{T} [GeV]; Entries",200,0,10);
784 TH1F* ptspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal p_{T} - MC; p_{T} [GeV]; Entries",200,0,10);
785 TH1F* ptspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background p_{T} - MC; p_{T} [GeV]; Entries",200,0,10);
786
787 ptspectrumMass->Sumw2();
788 ptspectrumSgn->Sumw2();
789 ptspectrumBkg->Sumw2();
790
791 ptspectrumMass->SetLineColor(6);
792 ptspectrumSgn->SetLineColor(2);
793 ptspectrumBkg->SetLineColor(4);
794
795 ptspectrumMass->SetMarkerStyle(20);
796 ptspectrumSgn->SetMarkerStyle(20);
797 ptspectrumBkg->SetMarkerStyle(20);
798 ptspectrumMass->SetMarkerSize(0.6);
799 ptspectrumSgn->SetMarkerSize(0.6);
800 ptspectrumBkg->SetMarkerSize(0.6);
801 ptspectrumMass->SetMarkerColor(6);
802 ptspectrumSgn->SetMarkerColor(2);
803 ptspectrumBkg->SetMarkerColor(4);
804
805 TH1F* ptallMass = (TH1F*)ptspectrumMass->Clone();
806 TH1F* ptallSgn = (TH1F*)ptspectrumSgn->Clone();
807 TH1F* ptallBkg = (TH1F*)ptspectrumBkg->Clone();
808
809 TH1F* ptpidMass = (TH1F*)ptspectrumMass->Clone();
810 TH1F* ptpidSgn = (TH1F*)ptspectrumSgn->Clone();
811 TH1F* ptpidBkg = (TH1F*)ptspectrumBkg->Clone();
812
813 fOutputAll->Add(ptallMass);
814 fOutputAll->Add(ptallSgn);
815 fOutputAll->Add(ptallBkg);
816
817 fOutputPID->Add(ptpidMass);
818 fOutputPID->Add(ptpidSgn);
819 fOutputPID->Add(ptpidBkg);
820
821 // eta spectra
822 nameMass="etaMass";
823 nameSgn="etaSgn";
824 nameBkg="etaBkg";
825
826 TH1F* etaspectrumMass = new TH1F(nameMass.Data(),"D^{*} #eta; #eta; Entries",200,-1,1);
827 TH1F* etaspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal #eta - MC; #eta; Entries",200,-1,1);
828 TH1F* etaspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background #eta - MC; #eta; Entries",200,-1,1);
829
830 etaspectrumMass->Sumw2();
831 etaspectrumSgn->Sumw2();
832 etaspectrumBkg->Sumw2();
833
834 etaspectrumMass->SetLineColor(6);
835 etaspectrumSgn->SetLineColor(2);
836 etaspectrumBkg->SetLineColor(4);
837
838 etaspectrumMass->SetMarkerStyle(20);
839 etaspectrumSgn->SetMarkerStyle(20);
840 etaspectrumBkg->SetMarkerStyle(20);
841 etaspectrumMass->SetMarkerSize(0.6);
842 etaspectrumSgn->SetMarkerSize(0.6);
843 etaspectrumBkg->SetMarkerSize(0.6);
844 etaspectrumMass->SetMarkerColor(6);
845 etaspectrumSgn->SetMarkerColor(2);
846 etaspectrumBkg->SetMarkerColor(4);
847
848 TH1F* etaallMass = (TH1F*)etaspectrumMass->Clone();
849 TH1F* etaallSgn = (TH1F*)etaspectrumSgn->Clone();
850 TH1F* etaallBkg = (TH1F*)etaspectrumBkg->Clone();
851
852 TH1F* etapidMass = (TH1F*)etaspectrumMass->Clone();
853 TH1F* etapidSgn = (TH1F*)etaspectrumSgn->Clone();
854 TH1F* etapidBkg = (TH1F*)etaspectrumBkg->Clone();
855
856 fOutputAll->Add(etaallMass);
857 fOutputAll->Add(etaallSgn);
858 fOutputAll->Add(etaallBkg);
859
860 fOutputPID->Add(etapidMass);
861 fOutputPID->Add(etapidSgn);
862 fOutputPID->Add(etapidBkg);
863
864 return;
865}
866//________________________________________________________________________
867void AliAnalysisTaskSEDStarSpectra::FillSpectrum(AliAODRecoCascadeHF *part, Int_t isDStar, AliRDHFCutsDStartoKpipi *cuts,Int_t isSel, TList *listout){
868 //
869 // Fill histos for D* spectrum
870 //
871
872 if(!isSel) return;
873
874 // D0 window
875 Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
876 Double_t invmassD0 = part->InvMassD0();
877 if (TMath::Abs(invmassD0-mPDGD0)>fD0Window) return;
878
879
880 Int_t ptbin=cuts->PtBin(part->Pt());
881 Double_t pt = part->Pt();
882 Double_t eta = part->Eta();
883
884 Double_t invmassDelta = part->DeltaInvMass();
885 Double_t invmassDstar = part->InvMassDstarKpipi();
886
887 TString fillthis="";
888 Bool_t massInRange=kFALSE;
889
890 Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
891
892 // delta M(Kpipi)-M(Kpi)
893 if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<fPeakWindow) massInRange=kTRUE;
894
895 if(fUseMCInfo) {
896 if(isDStar==1) {
897 fillthis="histD0Sgn_";
898 fillthis+=ptbin;
899 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
900 fillthis="histD0Sgn";
901 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
902 fillthis="histDstarSgn_";
903 fillthis+=ptbin;
904 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
905 fillthis="histDstarSgn";
906 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
907 fillthis="histDeltaSgn_";
908 fillthis+=ptbin;
909 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
910 fillthis="histDeltaSgn";
911 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
912 if (massInRange) {
913 fillthis="ptSgn";
914 ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
915 fillthis="etaSgn";
916 ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
917 }
918 }
919 else {//background
920 fillthis="histD0Bkg_";
921 fillthis+=ptbin;
922 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
923 fillthis="histD0Bkg";
924 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
925 fillthis="histDstarBkg_";
926 fillthis+=ptbin;
927 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
928 fillthis="histDstarBkg";
929 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
930 fillthis="histDeltaBkg_";
931 fillthis+=ptbin;
932 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
933 fillthis="histDeltaBkg";
934 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
935 if (massInRange) {
936 fillthis="ptBkg";
937 ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
938 fillthis="etaBkg";
939 ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
940 }
941 }
942 }
943 //no MC info, just cut selection
944 fillthis="histD0Mass_";
945 fillthis+=ptbin;
946 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
947 fillthis="histD0Mass";
948 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
949 fillthis="histDstarMass_";
950 fillthis+=ptbin;
951 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
952 fillthis="histDstarMass";
953 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
954 fillthis="histDeltaMass_";
955 fillthis+=ptbin;
956 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
957 fillthis="histDeltaMass";
958 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
959
960 if (massInRange) {
961 fillthis="ptMass";
962 ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
963 fillthis="etaMass";
964 ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
965 }
966
967 return;
968}
969//______________________________ side band background for D*___________________________________
970void AliAnalysisTaskSEDStarSpectra::SideBandBackground(AliAODRecoCascadeHF *part, AliRDHFCutsDStartoKpipi *cuts, Int_t isSel, TList *listout){
971
972 // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
973 // (expected detector resolution) on the left and right frm the D0 mass. Each band
974 // has a width of ~5 sigmas. Two band needed for opening angle considerations
975
976 if(!isSel) return;
977
978 Int_t ptbin=cuts->PtBin(part->Pt());
979
980 Bool_t massInRange=kFALSE;
981
982 // select the side bands intervall
983 Double_t invmassD0 = part->InvMassD0();
984 if(TMath::Abs(invmassD0-1.865)>4*fD0Window && TMath::Abs(invmassD0-1.865)<8*fD0Window){
985
986 // for pt and eta
987 Double_t invmassDelta = part->DeltaInvMass();
988 if (TMath::Abs(invmassDelta-0.14557)<fPeakWindow) massInRange=kTRUE;
989
990 TString fillthis="";
991 fillthis="histSideBandMass_";
992 fillthis+=ptbin;
993 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
994 fillthis="histSideBandMass";
995 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
996
997 }
998}
999//________________________________________________________________________________________________________________
1000void AliAnalysisTaskSEDStarSpectra::WrongSignForDStar(AliAODRecoCascadeHF *part, AliRDHFCutsDStartoKpipi *cuts, TList *listout){
1001 //
1002 // assign the wrong charge to the soft pion to create background
1003 //
1004 Int_t ptbin=cuts->PtBin(part->Pt());
1005
1006 Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1007 Double_t invmassD0 = part->InvMassD0();
1008 if (TMath::Abs(invmassD0-mPDGD0)>fD0Window) return;
1009
1010 AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)part->Get2Prong();
1011
1012 Int_t okD0WrongSign,okD0barWrongSign;
1013 Double_t wrongMassD0=0.;
1014
1015 Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kCandidate); //selected
1016 if (!isSelected){
1017 return;
1018 }
1019
1020 okD0WrongSign = 1;
1021 okD0barWrongSign = 1;
1022
1023 //if is D*+ than assume D0bar
1024 if(part->Charge()>0 && (isSelected ==1)) {
1025 okD0WrongSign = 0;
1026 }
1027 if(part->Charge()<0 && (isSelected ==2)){
1028 okD0barWrongSign = 0;
1029 }
1030
1031 // assign the wrong mass in case the cuts return both D0 and D0bar
1032 if(part->Charge()>0 && (isSelected ==3)) {
1033 okD0WrongSign = 0;
1034 } else if(part->Charge()<0 && (isSelected ==3)){
1035 okD0barWrongSign = 0;
1036 }
1037
1038 //wrong D0 inv mass
1039 if(okD0WrongSign!=0){
1040 wrongMassD0 = theD0particle->InvMassD0();
1041 }else if(okD0WrongSign==0){
1042 wrongMassD0 = theD0particle->InvMassD0bar();
1043 }
1044
1045 if(TMath::Abs(wrongMassD0-1.865)<fD0Window){
1046
1047 // wrong D* inv mass
1048 Double_t e[3];
1049 if (part->Charge()>0){
1050 e[0]=theD0particle->EProng(0,321);
1051 e[1]=theD0particle->EProng(1,211);
1052 }else{
1053 e[0]=theD0particle->EProng(0,211);
1054 e[1]=theD0particle->EProng(1,321);
1055 }
1056 e[2]=part->EProng(0,211);
1057
1058 Double_t esum = e[0]+e[1]+e[2];
1059 Double_t pds = part->P();
1060
1061 Double_t wrongMassDstar = TMath::Sqrt(esum*esum-pds*pds);
1062
1063 TString fillthis="";
1064 fillthis="histWrongSignMass_";
1065 fillthis+=ptbin;
1066 ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);
1067 fillthis="histWrongSignMass";
1068 ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);
1069
1070 }
1071}
1072//-------------------------------------------------------------------------------
1073Int_t AliAnalysisTaskSEDStarSpectra::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
1074 //
1075 // checking whether the mother of the particles come from a charm or a bottom quark
1076 //
1077
1078 Int_t pdgGranma = 0;
1079 Int_t mother = 0;
1080 mother = mcPartCandidate->GetMother();
1081 Int_t istep = 0;
1082 Int_t abspdgGranma =0;
1083 Bool_t isFromB=kFALSE;
1084 Bool_t isQuarkFound=kFALSE;
1085 while (mother >0 ){
1086 istep++;
1087 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1088 if (mcGranma){
1089 pdgGranma = mcGranma->GetPdgCode();
1090 abspdgGranma = TMath::Abs(pdgGranma);
1091 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1092 isFromB=kTRUE;
1093 }
1094 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1095 mother = mcGranma->GetMother();
1096 }else{
1097 AliError("Failed casting the mother particle!");
1098 break;
1099 }
1100 }
1101
1102 if(isFromB) return 5;
1103 else return 4;
1104}
1105//-------------------------------------------------------------------------------------
1106Float_t AliAnalysisTaskSEDStarSpectra::GetTrueImpactParameterD0(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partDp) const {
1107 // true impact parameter calculation
1108
1109 Double_t vtxTrue[3];
1110 mcHeader->GetVertex(vtxTrue);
1111 Double_t origD[3];
1112 partDp->XvYvZv(origD);
1113 Short_t charge=partDp->Charge();
1114 Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
1115 Int_t labelFirstDau = partDp->GetDaughter(0);
1116
1117 Int_t nDau=partDp->GetNDaughters();
1118
1119 Int_t theDau=0;
1120 if(nDau==2){
1121 for(Int_t iDau=0; iDau<2; iDau++){
1122 Int_t ind = labelFirstDau+iDau;
1123 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1124 if(!part){
1125 AliError("Daughter particle not found in MC array");
1126 return 99999.;
1127 }
1128 Int_t pdgCode=TMath::Abs(part->GetPdgCode());
1129 if(pdgCode==211 || pdgCode==321){
1130 pXdauTrue[theDau]=part->Px();
1131 pYdauTrue[theDau]=part->Py();
1132 pZdauTrue[theDau]=part->Pz();
1133 ++theDau;
1134 }
1135 }
1136 }
1137 if(theDau!=2){
1138 AliError("Wrong number of decay prongs");
1139 return 99999.;
1140 }
1141
1142 Double_t d0dummy[3]={0.,0.,0.};
1143 AliAODRecoDecayHF aodD0MC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
1144 return aodD0MC.ImpParXY();
1145
1146}
1147//______________________________________________________-
1148void AliAnalysisTaskSEDStarSpectra::CreateImpactParameterHistos(){
1149 // Histos for impact paramter study
1150
1151 Int_t nbins[3]={400,200,fNImpParBins};
1152 Double_t xmin[3]={1.75,0.,fLowerImpPar};
1153 Double_t xmax[3]={1.98,20.,fHigherImpPar};
1154
1155 fHistMassPtImpParTCDs[0]=new THnSparseF("hMassPtImpParAll",
1156 "Mass vs. pt vs.imppar - All",
1157 3,nbins,xmin,xmax);
1158 fHistMassPtImpParTCDs[1]=new THnSparseF("hMassPtImpParPrompt",
1159 "Mass vs. pt vs.imppar - promptD",
1160 3,nbins,xmin,xmax);
1161 fHistMassPtImpParTCDs[2]=new THnSparseF("hMassPtImpParBfeed",
1162 "Mass vs. pt vs.imppar - DfromB",
1163 3,nbins,xmin,xmax);
1164 fHistMassPtImpParTCDs[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1165 "Mass vs. pt vs.true imppar -DfromB",
1166 3,nbins,xmin,xmax);
1167 fHistMassPtImpParTCDs[4]=new THnSparseF("hMassPtImpParBkg",
1168 "Mass vs. pt vs.imppar - backgr.",
1169 3,nbins,xmin,xmax);
1170
1171 for(Int_t i=0; i<5;i++){
1172 fOutput->Add(fHistMassPtImpParTCDs[i]);
1173 }
1174}
1175