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