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