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