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