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