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