]> 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
Added Init() for AliNormalizationCounter
[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()));
38802708 554 fCounter->Init();
63147a75 555
556 if(fDoImpParDstar) CreateImpactParameterHistos();
557
558 PostData(1,fOutput);
559 PostData(2,fOutputAll);
560 PostData(3,fOutputPID);
561
562 return;
563}
564//___________________________________ hiostograms _______________________________________
565void AliAnalysisTaskSEDStarSpectra::DefineHistograms(){
566
567 fCEvents = new TH1F("fCEvents","conter",11,0,11);
568 fCEvents->SetStats(kTRUE);
569 fCEvents->GetXaxis()->SetTitle("1");
570 fCEvents->GetYaxis()->SetTitle("counts");
571 fOutput->Add(fCEvents);
572
573 fTrueDiff2 = new TH2F("DiffDstar_pt","True Reco diff vs pt",200,0,15,900,0,0.3);
574 fOutput->Add(fTrueDiff2);
575
576 fDeltaMassD1 = new TH1F("DeltaMassD1","delta mass d1",600,0,0.8);
577 fOutput->Add(fDeltaMassD1);
578
579 const Int_t nhist=14;
580 TString nameMass=" ", nameSgn=" ", nameBkg=" ";
581
582 for(Int_t i=-2;i<nhist;i++){
583 nameMass="histDeltaMass_";
584 nameMass+=i+1;
585 nameSgn="histDeltaSgn_";
586 nameSgn+=i+1;
587 nameBkg="histDeltaBkg_";
588 nameBkg+=i+1;
589
590 if (i==-2) {
591 nameMass="histDeltaMass";
592 nameSgn="histDeltaSgn";
593 nameBkg="histDeltaBkg";
594 }
595
596 TH1F* spectrumMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} invariant mass; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
597 TH1F* spectrumSgn = new TH1F(nameSgn.Data(), "D^{*}-D^{0} Signal invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
598 TH1F* spectrumBkg = new TH1F(nameBkg.Data(), "D^{*}-D^{0} Background invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
599
600 nameMass="histD0Mass_";
601 nameMass+=i+1;
602 nameSgn="histD0Sgn_";
603 nameSgn+=i+1;
604 nameBkg="histD0Bkg_";
605 nameBkg+=i+1;
606
607 if (i==-2) {
608 nameMass="histD0Mass";
609 nameSgn="histD0Sgn";
610 nameBkg="histD0Bkg";
611 }
612
613 TH1F* spectrumD0Mass = new TH1F(nameMass.Data(),"D^{0} invariant mass; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
614 TH1F* spectrumD0Sgn = new TH1F(nameSgn.Data(), "D^{0} Signal invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
615 TH1F* spectrumD0Bkg = new TH1F(nameBkg.Data(), "D^{0} Background invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
616
617 nameMass="histDstarMass_";
618 nameMass+=i+1;
619 nameSgn="histDstarSgn_";
620 nameSgn+=i+1;
621 nameBkg="histDstarBkg_";
622 nameBkg+=i+1;
623
624 if (i==-2) {
625 nameMass="histDstarMass";
626 nameSgn="histDstarSgn";
627 nameBkg="histDstarBkg";
628 }
629
630 TH1F* spectrumDstarMass = new TH1F(nameMass.Data(),"D^{*} invariant mass; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
631 TH1F* spectrumDstarSgn = new TH1F(nameSgn.Data(), "D^{*} Signal invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
632 TH1F* spectrumDstarBkg = new TH1F(nameBkg.Data(), "D^{*} Background invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
633
634 nameMass="histSideBandMass_";
635 nameMass+=i+1;
636 if (i==-2) {
637 nameMass="histSideBandMass";
638 }
639
640 TH1F* spectrumSideBandMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} sideband mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
641
642 nameMass="histWrongSignMass_";
643 nameMass+=i+1;
644 if (i==-2) {
645 nameMass="histWrongSignMass";
646 }
647
648 TH1F* spectrumWrongSignMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} wrongsign mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
649
650
651 spectrumMass->Sumw2();
652 spectrumSgn->Sumw2();
653 spectrumBkg->Sumw2();
654
655 spectrumMass->SetLineColor(6);
656 spectrumSgn->SetLineColor(2);
657 spectrumBkg->SetLineColor(4);
658
659 spectrumMass->SetMarkerStyle(20);
660 spectrumSgn->SetMarkerStyle(20);
661 spectrumBkg->SetMarkerStyle(20);
662 spectrumMass->SetMarkerSize(0.6);
663 spectrumSgn->SetMarkerSize(0.6);
664 spectrumBkg->SetMarkerSize(0.6);
665 spectrumMass->SetMarkerColor(6);
666 spectrumSgn->SetMarkerColor(2);
667 spectrumBkg->SetMarkerColor(4);
668
669 spectrumD0Mass->Sumw2();
670 spectrumD0Sgn->Sumw2();
671 spectrumD0Bkg->Sumw2();
672
673 spectrumD0Mass->SetLineColor(6);
674 spectrumD0Sgn->SetLineColor(2);
675 spectrumD0Bkg->SetLineColor(4);
676
677 spectrumD0Mass->SetMarkerStyle(20);
678 spectrumD0Sgn->SetMarkerStyle(20);
679 spectrumD0Bkg->SetMarkerStyle(20);
680 spectrumD0Mass->SetMarkerSize(0.6);
681 spectrumD0Sgn->SetMarkerSize(0.6);
682 spectrumD0Bkg->SetMarkerSize(0.6);
683 spectrumD0Mass->SetMarkerColor(6);
684 spectrumD0Sgn->SetMarkerColor(2);
685 spectrumD0Bkg->SetMarkerColor(4);
686
687 spectrumDstarMass->Sumw2();
688 spectrumDstarSgn->Sumw2();
689 spectrumDstarBkg->Sumw2();
690
691 spectrumDstarMass->SetLineColor(6);
692 spectrumDstarSgn->SetLineColor(2);
693 spectrumDstarBkg->SetLineColor(4);
694
695 spectrumDstarMass->SetMarkerStyle(20);
696 spectrumDstarSgn->SetMarkerStyle(20);
697 spectrumDstarBkg->SetMarkerStyle(20);
698 spectrumDstarMass->SetMarkerSize(0.6);
699 spectrumDstarSgn->SetMarkerSize(0.6);
700 spectrumDstarBkg->SetMarkerSize(0.6);
701 spectrumDstarMass->SetMarkerColor(6);
702 spectrumDstarSgn->SetMarkerColor(2);
703 spectrumDstarBkg->SetMarkerColor(4);
704
705 spectrumSideBandMass->Sumw2();
706 spectrumSideBandMass->SetLineColor(4);
707 spectrumSideBandMass->SetMarkerStyle(20);
708 spectrumSideBandMass->SetMarkerSize(0.6);
709 spectrumSideBandMass->SetMarkerColor(4);
710
711 spectrumWrongSignMass->Sumw2();
712 spectrumWrongSignMass->SetLineColor(4);
713 spectrumWrongSignMass->SetMarkerStyle(20);
714 spectrumWrongSignMass->SetMarkerSize(0.6);
715 spectrumWrongSignMass->SetMarkerColor(4);
716
717 TH1F* allMass = (TH1F*)spectrumMass->Clone();
718 TH1F* allSgn = (TH1F*)spectrumSgn->Clone();
719 TH1F* allBkg = (TH1F*)spectrumBkg->Clone();
720
721 TH1F* pidMass = (TH1F*)spectrumMass->Clone();
722 TH1F* pidSgn = (TH1F*)spectrumSgn->Clone();
723 TH1F* pidBkg = (TH1F*)spectrumBkg->Clone();
724
725 fOutputAll->Add(allMass);
726 fOutputAll->Add(allSgn);
727 fOutputAll->Add(allBkg);
728
729 fOutputPID->Add(pidMass);
730 fOutputPID->Add(pidSgn);
731 fOutputPID->Add(pidBkg);
732
733 TH1F* allD0Mass = (TH1F*)spectrumD0Mass->Clone();
734 TH1F* allD0Sgn = (TH1F*)spectrumD0Sgn->Clone();
735 TH1F* allD0Bkg = (TH1F*)spectrumD0Bkg->Clone();
736
737 TH1F* pidD0Mass = (TH1F*)spectrumD0Mass->Clone();
738 TH1F* pidD0Sgn = (TH1F*)spectrumD0Sgn->Clone();
739 TH1F* pidD0Bkg = (TH1F*)spectrumD0Bkg->Clone();
740
741 fOutputAll->Add(allD0Mass);
742 fOutputAll->Add(allD0Sgn);
743 fOutputAll->Add(allD0Bkg);
744
745 fOutputPID->Add(pidD0Mass);
746 fOutputPID->Add(pidD0Sgn);
747 fOutputPID->Add(pidD0Bkg);
748
749 TH1F* allDstarMass = (TH1F*)spectrumDstarMass->Clone();
750 TH1F* allDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
751 TH1F* allDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
752
753 TH1F* pidDstarMass = (TH1F*)spectrumDstarMass->Clone();
754 TH1F* pidDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
755 TH1F* pidDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
756
757 fOutputAll->Add(allDstarMass);
758 fOutputAll->Add(allDstarSgn);
759 fOutputAll->Add(allDstarBkg);
760
761 fOutputPID->Add(pidDstarMass);
762 fOutputPID->Add(pidDstarSgn);
763 fOutputPID->Add(pidDstarBkg);
764
765 TH1F* allSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
766 TH1F* pidSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
767
768 fOutputAll->Add(allSideBandMass);
769 fOutputPID->Add(pidSideBandMass);
770
771 TH1F* allWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
772 TH1F* pidWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
773
774 fOutputAll->Add(allWrongSignMass);
775 fOutputPID->Add(pidWrongSignMass);
776
777 }
778
779 // pt spectra
780 nameMass="ptMass";
781 nameSgn="ptSgn";
782 nameBkg="ptBkg";
783
784 TH1F* ptspectrumMass = new TH1F(nameMass.Data(),"D^{*} p_{T}; p_{T} [GeV]; Entries",200,0,10);
785 TH1F* ptspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal p_{T} - MC; p_{T} [GeV]; Entries",200,0,10);
786 TH1F* ptspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background p_{T} - MC; p_{T} [GeV]; Entries",200,0,10);
787
788 ptspectrumMass->Sumw2();
789 ptspectrumSgn->Sumw2();
790 ptspectrumBkg->Sumw2();
791
792 ptspectrumMass->SetLineColor(6);
793 ptspectrumSgn->SetLineColor(2);
794 ptspectrumBkg->SetLineColor(4);
795
796 ptspectrumMass->SetMarkerStyle(20);
797 ptspectrumSgn->SetMarkerStyle(20);
798 ptspectrumBkg->SetMarkerStyle(20);
799 ptspectrumMass->SetMarkerSize(0.6);
800 ptspectrumSgn->SetMarkerSize(0.6);
801 ptspectrumBkg->SetMarkerSize(0.6);
802 ptspectrumMass->SetMarkerColor(6);
803 ptspectrumSgn->SetMarkerColor(2);
804 ptspectrumBkg->SetMarkerColor(4);
805
806 TH1F* ptallMass = (TH1F*)ptspectrumMass->Clone();
807 TH1F* ptallSgn = (TH1F*)ptspectrumSgn->Clone();
808 TH1F* ptallBkg = (TH1F*)ptspectrumBkg->Clone();
809
810 TH1F* ptpidMass = (TH1F*)ptspectrumMass->Clone();
811 TH1F* ptpidSgn = (TH1F*)ptspectrumSgn->Clone();
812 TH1F* ptpidBkg = (TH1F*)ptspectrumBkg->Clone();
813
814 fOutputAll->Add(ptallMass);
815 fOutputAll->Add(ptallSgn);
816 fOutputAll->Add(ptallBkg);
817
818 fOutputPID->Add(ptpidMass);
819 fOutputPID->Add(ptpidSgn);
820 fOutputPID->Add(ptpidBkg);
821
822 // eta spectra
823 nameMass="etaMass";
824 nameSgn="etaSgn";
825 nameBkg="etaBkg";
826
827 TH1F* etaspectrumMass = new TH1F(nameMass.Data(),"D^{*} #eta; #eta; Entries",200,-1,1);
828 TH1F* etaspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal #eta - MC; #eta; Entries",200,-1,1);
829 TH1F* etaspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background #eta - MC; #eta; Entries",200,-1,1);
830
831 etaspectrumMass->Sumw2();
832 etaspectrumSgn->Sumw2();
833 etaspectrumBkg->Sumw2();
834
835 etaspectrumMass->SetLineColor(6);
836 etaspectrumSgn->SetLineColor(2);
837 etaspectrumBkg->SetLineColor(4);
838
839 etaspectrumMass->SetMarkerStyle(20);
840 etaspectrumSgn->SetMarkerStyle(20);
841 etaspectrumBkg->SetMarkerStyle(20);
842 etaspectrumMass->SetMarkerSize(0.6);
843 etaspectrumSgn->SetMarkerSize(0.6);
844 etaspectrumBkg->SetMarkerSize(0.6);
845 etaspectrumMass->SetMarkerColor(6);
846 etaspectrumSgn->SetMarkerColor(2);
847 etaspectrumBkg->SetMarkerColor(4);
848
849 TH1F* etaallMass = (TH1F*)etaspectrumMass->Clone();
850 TH1F* etaallSgn = (TH1F*)etaspectrumSgn->Clone();
851 TH1F* etaallBkg = (TH1F*)etaspectrumBkg->Clone();
852
853 TH1F* etapidMass = (TH1F*)etaspectrumMass->Clone();
854 TH1F* etapidSgn = (TH1F*)etaspectrumSgn->Clone();
855 TH1F* etapidBkg = (TH1F*)etaspectrumBkg->Clone();
856
857 fOutputAll->Add(etaallMass);
858 fOutputAll->Add(etaallSgn);
859 fOutputAll->Add(etaallBkg);
860
861 fOutputPID->Add(etapidMass);
862 fOutputPID->Add(etapidSgn);
863 fOutputPID->Add(etapidBkg);
864
865 return;
866}
867//________________________________________________________________________
868void AliAnalysisTaskSEDStarSpectra::FillSpectrum(AliAODRecoCascadeHF *part, Int_t isDStar, AliRDHFCutsDStartoKpipi *cuts,Int_t isSel, TList *listout){
869 //
870 // Fill histos for D* spectrum
871 //
872
873 if(!isSel) return;
874
875 // D0 window
876 Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
877 Double_t invmassD0 = part->InvMassD0();
878 if (TMath::Abs(invmassD0-mPDGD0)>fD0Window) return;
879
880
881 Int_t ptbin=cuts->PtBin(part->Pt());
882 Double_t pt = part->Pt();
883 Double_t eta = part->Eta();
884
885 Double_t invmassDelta = part->DeltaInvMass();
886 Double_t invmassDstar = part->InvMassDstarKpipi();
887
888 TString fillthis="";
889 Bool_t massInRange=kFALSE;
890
891 Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
892
893 // delta M(Kpipi)-M(Kpi)
894 if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<fPeakWindow) massInRange=kTRUE;
895
896 if(fUseMCInfo) {
897 if(isDStar==1) {
898 fillthis="histD0Sgn_";
899 fillthis+=ptbin;
900 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
901 fillthis="histD0Sgn";
902 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
903 fillthis="histDstarSgn_";
904 fillthis+=ptbin;
905 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
906 fillthis="histDstarSgn";
907 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
908 fillthis="histDeltaSgn_";
909 fillthis+=ptbin;
910 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
911 fillthis="histDeltaSgn";
912 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
913 if (massInRange) {
914 fillthis="ptSgn";
915 ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
916 fillthis="etaSgn";
917 ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
918 }
919 }
920 else {//background
921 fillthis="histD0Bkg_";
922 fillthis+=ptbin;
923 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
924 fillthis="histD0Bkg";
925 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
926 fillthis="histDstarBkg_";
927 fillthis+=ptbin;
928 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
929 fillthis="histDstarBkg";
930 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
931 fillthis="histDeltaBkg_";
932 fillthis+=ptbin;
933 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
934 fillthis="histDeltaBkg";
935 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
936 if (massInRange) {
937 fillthis="ptBkg";
938 ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
939 fillthis="etaBkg";
940 ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
941 }
942 }
943 }
944 //no MC info, just cut selection
945 fillthis="histD0Mass_";
946 fillthis+=ptbin;
947 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
948 fillthis="histD0Mass";
949 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
950 fillthis="histDstarMass_";
951 fillthis+=ptbin;
952 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
953 fillthis="histDstarMass";
954 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
955 fillthis="histDeltaMass_";
956 fillthis+=ptbin;
957 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
958 fillthis="histDeltaMass";
959 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
960
961 if (massInRange) {
962 fillthis="ptMass";
963 ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
964 fillthis="etaMass";
965 ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
966 }
967
968 return;
969}
970//______________________________ side band background for D*___________________________________
971void AliAnalysisTaskSEDStarSpectra::SideBandBackground(AliAODRecoCascadeHF *part, AliRDHFCutsDStartoKpipi *cuts, Int_t isSel, TList *listout){
972
973 // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
974 // (expected detector resolution) on the left and right frm the D0 mass. Each band
975 // has a width of ~5 sigmas. Two band needed for opening angle considerations
976
977 if(!isSel) return;
978
979 Int_t ptbin=cuts->PtBin(part->Pt());
980
981 Bool_t massInRange=kFALSE;
982
983 // select the side bands intervall
984 Double_t invmassD0 = part->InvMassD0();
985 if(TMath::Abs(invmassD0-1.865)>4*fD0Window && TMath::Abs(invmassD0-1.865)<8*fD0Window){
986
987 // for pt and eta
988 Double_t invmassDelta = part->DeltaInvMass();
989 if (TMath::Abs(invmassDelta-0.14557)<fPeakWindow) massInRange=kTRUE;
990
991 TString fillthis="";
992 fillthis="histSideBandMass_";
993 fillthis+=ptbin;
994 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
995 fillthis="histSideBandMass";
996 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
997
998 }
999}
1000//________________________________________________________________________________________________________________
1001void AliAnalysisTaskSEDStarSpectra::WrongSignForDStar(AliAODRecoCascadeHF *part, AliRDHFCutsDStartoKpipi *cuts, TList *listout){
1002 //
1003 // assign the wrong charge to the soft pion to create background
1004 //
1005 Int_t ptbin=cuts->PtBin(part->Pt());
1006
1007 Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1008 Double_t invmassD0 = part->InvMassD0();
1009 if (TMath::Abs(invmassD0-mPDGD0)>fD0Window) return;
1010
1011 AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)part->Get2Prong();
1012
1013 Int_t okD0WrongSign,okD0barWrongSign;
1014 Double_t wrongMassD0=0.;
1015
1016 Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kCandidate); //selected
1017 if (!isSelected){
1018 return;
1019 }
1020
1021 okD0WrongSign = 1;
1022 okD0barWrongSign = 1;
1023
1024 //if is D*+ than assume D0bar
1025 if(part->Charge()>0 && (isSelected ==1)) {
1026 okD0WrongSign = 0;
1027 }
1028 if(part->Charge()<0 && (isSelected ==2)){
1029 okD0barWrongSign = 0;
1030 }
1031
1032 // assign the wrong mass in case the cuts return both D0 and D0bar
1033 if(part->Charge()>0 && (isSelected ==3)) {
1034 okD0WrongSign = 0;
1035 } else if(part->Charge()<0 && (isSelected ==3)){
1036 okD0barWrongSign = 0;
1037 }
1038
1039 //wrong D0 inv mass
1040 if(okD0WrongSign!=0){
1041 wrongMassD0 = theD0particle->InvMassD0();
1042 }else if(okD0WrongSign==0){
1043 wrongMassD0 = theD0particle->InvMassD0bar();
1044 }
1045
1046 if(TMath::Abs(wrongMassD0-1.865)<fD0Window){
1047
1048 // wrong D* inv mass
1049 Double_t e[3];
1050 if (part->Charge()>0){
1051 e[0]=theD0particle->EProng(0,321);
1052 e[1]=theD0particle->EProng(1,211);
1053 }else{
1054 e[0]=theD0particle->EProng(0,211);
1055 e[1]=theD0particle->EProng(1,321);
1056 }
1057 e[2]=part->EProng(0,211);
1058
1059 Double_t esum = e[0]+e[1]+e[2];
1060 Double_t pds = part->P();
1061
1062 Double_t wrongMassDstar = TMath::Sqrt(esum*esum-pds*pds);
1063
1064 TString fillthis="";
1065 fillthis="histWrongSignMass_";
1066 fillthis+=ptbin;
1067 ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);
1068 fillthis="histWrongSignMass";
1069 ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);
1070
1071 }
1072}
1073//-------------------------------------------------------------------------------
1074Int_t AliAnalysisTaskSEDStarSpectra::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
1075 //
1076 // checking whether the mother of the particles come from a charm or a bottom quark
1077 //
1078
1079 Int_t pdgGranma = 0;
1080 Int_t mother = 0;
1081 mother = mcPartCandidate->GetMother();
1082 Int_t istep = 0;
1083 Int_t abspdgGranma =0;
1084 Bool_t isFromB=kFALSE;
1085 Bool_t isQuarkFound=kFALSE;
1086 while (mother >0 ){
1087 istep++;
1088 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1089 if (mcGranma){
1090 pdgGranma = mcGranma->GetPdgCode();
1091 abspdgGranma = TMath::Abs(pdgGranma);
1092 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1093 isFromB=kTRUE;
1094 }
1095 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1096 mother = mcGranma->GetMother();
1097 }else{
1098 AliError("Failed casting the mother particle!");
1099 break;
1100 }
1101 }
1102
1103 if(isFromB) return 5;
1104 else return 4;
1105}
1106//-------------------------------------------------------------------------------------
1107Float_t AliAnalysisTaskSEDStarSpectra::GetTrueImpactParameterD0(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partDp) const {
1108 // true impact parameter calculation
1109
1110 Double_t vtxTrue[3];
1111 mcHeader->GetVertex(vtxTrue);
1112 Double_t origD[3];
1113 partDp->XvYvZv(origD);
1114 Short_t charge=partDp->Charge();
1115 Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
1116 Int_t labelFirstDau = partDp->GetDaughter(0);
1117
1118 Int_t nDau=partDp->GetNDaughters();
1119
1120 Int_t theDau=0;
1121 if(nDau==2){
1122 for(Int_t iDau=0; iDau<2; iDau++){
1123 Int_t ind = labelFirstDau+iDau;
1124 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1125 if(!part){
1126 AliError("Daughter particle not found in MC array");
1127 return 99999.;
1128 }
1129 Int_t pdgCode=TMath::Abs(part->GetPdgCode());
1130 if(pdgCode==211 || pdgCode==321){
1131 pXdauTrue[theDau]=part->Px();
1132 pYdauTrue[theDau]=part->Py();
1133 pZdauTrue[theDau]=part->Pz();
1134 ++theDau;
1135 }
1136 }
1137 }
1138 if(theDau!=2){
1139 AliError("Wrong number of decay prongs");
1140 return 99999.;
1141 }
1142
1143 Double_t d0dummy[3]={0.,0.,0.};
1144 AliAODRecoDecayHF aodD0MC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
1145 return aodD0MC.ImpParXY();
1146
1147}
1148//______________________________________________________-
1149void AliAnalysisTaskSEDStarSpectra::CreateImpactParameterHistos(){
1150 // Histos for impact paramter study
1151
1152 Int_t nbins[3]={400,200,fNImpParBins};
1153 Double_t xmin[3]={1.75,0.,fLowerImpPar};
1154 Double_t xmax[3]={1.98,20.,fHigherImpPar};
1155
1156 fHistMassPtImpParTCDs[0]=new THnSparseF("hMassPtImpParAll",
1157 "Mass vs. pt vs.imppar - All",
1158 3,nbins,xmin,xmax);
1159 fHistMassPtImpParTCDs[1]=new THnSparseF("hMassPtImpParPrompt",
1160 "Mass vs. pt vs.imppar - promptD",
1161 3,nbins,xmin,xmax);
1162 fHistMassPtImpParTCDs[2]=new THnSparseF("hMassPtImpParBfeed",
1163 "Mass vs. pt vs.imppar - DfromB",
1164 3,nbins,xmin,xmax);
1165 fHistMassPtImpParTCDs[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1166 "Mass vs. pt vs.true imppar -DfromB",
1167 3,nbins,xmin,xmax);
1168 fHistMassPtImpParTCDs[4]=new THnSparseF("hMassPtImpParBkg",
1169 "Mass vs. pt vs.imppar - backgr.",
1170 3,nbins,xmin,xmax);
1171
1172 for(Int_t i=0; i<5;i++){
1173 fOutput->Add(fHistMassPtImpParTCDs[i]);
1174 }
1175}
1176