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