[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDStarSpectra.cxx

/**************************************************************************\r * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r *                                                                        *\r * Author: The ALICE Off-line Project.                                    *\r * Contributors are mentioned in the code where appropriate.              *\r *                                                                        *\r * Permission to use, copy, modify and distribute this software and its   *\r * documentation strictly for non-commercial purposes is hereby granted   *\r * without fee, provided that the above copyright notice appears in all   *\r * copies and that both the copyright notice and this permission notice   *\r * appeuear in the supporting documentation. The authors make no claims     *\r * about the suitability of this software for any purpose. It is          *\r * provided "as is" without express or implied warranty.                  *\r **************************************************************************/\r\r/* $Id$ */\r\r//\r//\r//                  Base class for DStar Analysis\r//\r//\r//  The D* spectra study is done in pt bins:\r//  [0,0.5] [0.5,1] [1,2] [2,3] [3,4] [4,5] [5,6] [6,7] [7,8],\r//  [8,10],[10,12], [12,16], [16,20] and [20,24]\r//\r//  Cuts arew centralized in AliRDHFCutsDStartoKpipi\r//  Side Band and like sign background are implemented in the macro\r//\r//-----------------------------------------------------------------------\r//\r//                         Author A.Grelli \r//              ERC-QGP Utrecht University - a.grelli@uu.nl,\r//                         Author Y.Wang\r//        University of Heidelberg - yifei@physi.uni-heidelberg.de\r//                         Author C.Ivan \r//             ERC-QGP Utrecht University - c.ivan@uu.nl,\r//\r//-----------------------------------------------------------------------\r\r#include <TSystem.h>\r#include <TParticle.h>\r#include <TH1I.h>\r#include "TROOT.h"\r#include <TDatabasePDG.h>\r#include <AliAnalysisDataSlot.h>\r#include <AliAnalysisDataContainer.h>\r#include "AliRDHFCutsDStartoKpipi.h"\r#include "AliStack.h"\r#include "AliMCEvent.h"\r#include "AliAnalysisManager.h"\r#include "AliAODMCHeader.h"\r#include "AliAODHandler.h"\r#include "AliLog.h"\r#include "AliAODVertex.h"\r#include "AliAODRecoDecay.h"\r#include "AliAODRecoDecayHF.h"\r#include "AliAODRecoCascadeHF.h"\r#include "AliAODRecoDecayHF2Prong.h"\r#include "AliAnalysisVertexingHF.h"\r#include "AliESDtrack.h"\r#include "AliAODMCParticle.h"\r#include "AliAnalysisTaskSE.h"\r#include "AliAnalysisTaskSEDStarSpectra.h"\r#include "AliNormalizationCounter.h"\r\rClassImp(AliAnalysisTaskSEDStarSpectra)\r\r//__________________________________________________________________________\rAliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra():  \r  AliAnalysisTaskSE(),\r  fEvents(0),\r  fAnalysis(0),\r  fD0Window(0),\r  fPeakWindow(0),\r  fUseMCInfo(kFALSE),\r  fDoSearch(kFALSE),\r  fOutput(0),\r  fOutputAll(0),\r  fOutputPID(0),\r  fNSigma(3),\r  fCuts(0),\r  fCEvents(0),     \r  fTrueDiff2(0),\r  fDeltaMassD1(0),\r  fCounter(0),\r  fDoImpParDstar(kFALSE),\r  fNImpParBins(400),\r  fLowerImpPar(-2000.),\r  fHigherImpPar(2000.)\r{\r  //\r  // Default ctor\r  //\r}\r//___________________________________________________________________________\rAliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra(const Char_t* name, AliRDHFCutsDStartoKpipi* cuts) :\r  AliAnalysisTaskSE(name),\r  fEvents(0),\r  fAnalysis(0),\r  fD0Window(0),\r  fPeakWindow(0),\r  fUseMCInfo(kFALSE),\r  fDoSearch(kFALSE),\r  fOutput(0),\r  fOutputAll(0),\r  fOutputPID(0),\r  fNSigma(3),\r  fCuts(0),\r  fCEvents(0),     \r  fTrueDiff2(0),\r  fDeltaMassD1(0),\r  fCounter(0),\r  fDoImpParDstar(kFALSE),\r  fNImpParBins(400),\r  fLowerImpPar(-2000.),\r  fHigherImpPar(2000.)\r{\r  //\r  // Constructor. Initialization of Inputs and Outputs\r  //\r  Info("AliAnalysisTaskSEDStarSpectra","Calling Constructor");\r\r  fCuts=cuts;\r\r  DefineOutput(1,TList::Class());  //conters\r  DefineOutput(2,TList::Class());  //All Entries output\r  DefineOutput(3,TList::Class());  //3sigma PID output\r  DefineOutput(4,AliRDHFCutsDStartoKpipi::Class());   //My private output\r  DefineOutput(5,AliNormalizationCounter::Class());   // normalization\r}\r\r//___________________________________________________________________________\rAliAnalysisTaskSEDStarSpectra::~AliAnalysisTaskSEDStarSpectra() {\r  //\r  // destructor\r  //\r  Info("~AliAnalysisTaskSEDStarSpectra","Calling Destructor");\r  \r  if (fOutput) {\r    delete fOutput;\r    fOutput = 0;\r  }\r  if (fOutputAll) {\r    delete fOutputAll;\r    fOutputAll = 0;\r  }\r  if (fOutputPID) {\r    delete fOutputPID;\r    fOutputPID = 0;\r  }\r  if (fCuts) {\r    delete fCuts;\r    fCuts = 0;\r  }\r  if(fCEvents){\r    delete fCEvents;\r    fCEvents =0;\r  }\r  if(fDeltaMassD1){\r    delete fDeltaMassD1;\r    fDeltaMassD1 =0;\r  }\r  for(Int_t i=0; i<5; i++){\r    delete fHistMassPtImpParTCDs[i];\r  }\r}\r//_________________________________________________\rvoid AliAnalysisTaskSEDStarSpectra::Init(){\r  //\r  // Initialization\r  //\r\r  if(fDebug > 1) printf("AnalysisTaskSEDStarSpectra::Init() \n");\r   AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);\r  // Post the data\r  PostData(4,copyfCuts);\r\r  return;\r}\r\r//_________________________________________________\rvoid AliAnalysisTaskSEDStarSpectra::UserExec(Option_t *)\r{\r  // user exec\r  if (!fInputEvent) {\r    Error("UserExec","NO EVENT FOUND!");\r    return;\r  }\r\r  fEvents++;\r\r  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);\r  TClonesArray *arrayDStartoD0pi=0;\r\r  fCEvents->Fill(1);\r\r  if(!aodEvent && AODEvent() && IsStandardAOD()) {\r    // In case there is an AOD handler writing a standard AOD, use the AOD \r    // event in memory rather than the input (ESD) event.    \r    aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());\r    // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)\r    // have to taken from the AOD event hold by the AliAODExtension\r    AliAODHandler* aodHandler = (AliAODHandler*) \r      ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());\r    if(aodHandler->GetExtensions()) {\r      AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");\r      AliAODEvent *aodFromExt = ext->GetAOD();\r      arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");\r    }\r  } else {\r    arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");\r  }\r\r  // fix for temporary bug in ESDfilter \r  // the AODs with null vertex pointer didn't pass the PhysSel\r  if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;\r  fCEvents->Fill(2);\r\r  fCounter->StoreEvent(aodEvent,fCuts,fUseMCInfo);\r\r  // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD\r  TString trigclass=aodEvent->GetFiredTriggerClasses();\r  if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")||trigclass.Contains("C0SMH-B-NOPF-ALL")) fCEvents->Fill(5);\r\r  if(!fCuts->IsEventSelected(aodEvent)) {\r    if(fCuts->GetWhyRejection()==6) // rejected for Z vertex\r      fCEvents->Fill(6);\r    return;\r  }\r\r  Bool_t isEvSel=fCuts->IsEventSelected(aodEvent);\r  if(!isEvSel) return;\r\r  // Load the event\r  AliInfo(Form("Event %d",fEvents));\r  if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));\r\r  // counters for efficiencies\r  Int_t icountReco = 0;\r  \r  //D* and D0 prongs needed to MatchToMC method\r  Int_t pdgDgDStartoD0pi[2]={421,211};\r  Int_t pdgDgD0toKpi[2]={321,211};\r\r  // AOD primary vertex\r  AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();\r  if(!vtx1) return;\r  if(vtx1->GetNContributors()<1) return;\r\r  if (!arrayDStartoD0pi){\r    AliInfo("Could not find array of HF vertices, skipping the event");\r    return;\r  }else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast())); \r\r  Int_t nSelectedAna =0;\r  Int_t nSelectedProd =0;\r\r  // loop over the tracks to search for candidates soft pion \r  for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {\r\r    // D* candidates and D0 from D*\r    AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);\r    if(!dstarD0pi->GetSecondaryVtx()) continue;\r    AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();\r    if (!theD0particle) continue;\r    \r    Int_t isDStar = 0;   \r    TClonesArray *mcArray = 0;\r    AliAODMCHeader *mcHeader=0;\r\r    Bool_t isPrimary=kTRUE;\r    Float_t truePt=0.;\r    Float_t xDecay=0.;\r    Float_t yDecay=0.;\r    Float_t zDecay=0.;\r    Float_t pdgCode=-2;\r    Float_t trueImpParXY=0.;\r\r    // mc analysis \r    if(fUseMCInfo){\r    //MC array need for maching\r      mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));\r      if (!mcArray) {\r	AliError("Could not find Monte-Carlo in AOD");\r	return;\r      }\r      // load MC header\r      mcHeader =  (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());\r      if(!mcHeader) {\r	printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");\r	return;\r      }\r      // find associated MC particle for D* ->D0toKpi\r      Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);\r      if(mcLabel>=0){\r\r	AliAODMCParticle *partDSt = (AliAODMCParticle*)mcArray->At(mcLabel);\r	Int_t checkOrigin = CheckOrigin(mcArray,partDSt);\r       	if(checkOrigin==5) isPrimary=kFALSE;\r	AliAODMCParticle *dg0 = (AliAODMCParticle*)mcArray->At(partDSt->GetDaughter(0));\r	AliAODMCParticle *dg0_1 = (AliAODMCParticle*)mcArray->At(dg0->GetDaughter(0));\r	truePt=dg0->Pt();\r	xDecay=dg0_1->Xv();	  \r	yDecay=dg0_1->Yv();	  \r	zDecay=dg0_1->Zv();\r	pdgCode=TMath::Abs(partDSt->GetPdgCode());\r	if(!isPrimary){\r	  trueImpParXY=GetTrueImpactParameterD0(mcHeader,mcArray,dg0)*1000.;\r	}\r	isDStar = 1;\r      }else{\r	pdgCode=-1;\r      }\r    }\r    \r    Int_t ptbin=fCuts->PtBin(dstarD0pi->Pt());\r    \r    // quality selction on tracks and region of interest\r    Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks\r    if(!isTkSelected) continue;\r\r    if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;    \r\r\r    //histos for impact par studies - D0!!!\r    Double_t ptCand = dstarD0pi->Get2Prong()->Pt();\r    Double_t invMass=dstarD0pi->InvMassD0();\r    Double_t impparXY=dstarD0pi->Get2Prong()->ImpParXY()*10000.;\r\r    Double_t arrayForSparse[3]={invMass,ptCand,impparXY};\r    Double_t arrayForSparseTrue[3]={invMass,ptCand,trueImpParXY};\r   \r  // set the D0 search window bin by bin - useful to calculate side band bkg\r    if (ptbin==0){\r      if(fAnalysis==1){\r	fD0Window=0.035;\r	fPeakWindow=0.03;\r      }else{\r	fD0Window=0.020;\r	fPeakWindow=0.0018;\r      }\r    }\r    if (ptbin==1){\r      if(fAnalysis==1){\r	fD0Window=0.035;\r	fPeakWindow=0.03;\r      }else{\r	fD0Window=0.020;\r	fPeakWindow=0.0018;\r      }\r    }\r    if (ptbin==2){\r      if(fAnalysis==1){\r	fD0Window=0.035;\r	fPeakWindow=0.03;\r      }else{\r	fD0Window=0.020;\r	fPeakWindow=0.0018;\r      }\r    }\r    if (ptbin==3){\r      if(fAnalysis==1){\r	fD0Window=0.035;\r	fPeakWindow=0.03;\r      }else{\r	fD0Window=0.022;\r	fPeakWindow=0.0016;\r      }\r    }\r    if (ptbin==4){ \r      if(fAnalysis==1){\r	fD0Window=0.035;\r	fPeakWindow=0.03;\r      }else{\r	fD0Window=0.026;\r	fPeakWindow=0.0014;\r      }\r    }\r    if (ptbin==5){\r      if(fAnalysis==1){\r	fD0Window=0.045;\r	fPeakWindow=0.03;\r      }else{\r	fD0Window=0.026;\r	fPeakWindow=0.0014;\r      }\r    } \r   if (ptbin==6){\r      if(fAnalysis==1){\r	fD0Window=0.045;\r	fPeakWindow=0.03;\r      }else{\r	fD0Window=0.026;\r	fPeakWindow=0.006;\r      }\r    } \r    if (ptbin==7){\r      if(fAnalysis==1){\r	fD0Window=0.055;\r	fPeakWindow=0.03;\r      }else{\r	fD0Window=0.026;\r	fPeakWindow=0.006;\r      }\r    }\r    if (ptbin>7){\r      if(fAnalysis==1){\r	fD0Window=0.074;\r	fPeakWindow=0.03;\r      }else{\r	fD0Window=0.026;\r	fPeakWindow=0.006;\r      }\r    }\r    \r    nSelectedProd++;\r    nSelectedAna++;\r    \r    // check that we are close to signal in the DeltaM - here to save time for PbPb\r    Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();  \r    Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();\r    Double_t invmassDelta = dstarD0pi->DeltaInvMass();\r    \r    if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))>fPeakWindow) continue;\r    \r    \r    Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected\r    \r    // after cuts\r    if(fDoImpParDstar && isSelected){\r      fHistMassPtImpParTCDs[0]->Fill(arrayForSparse);\r      if(isPrimary) fHistMassPtImpParTCDs[1]->Fill(arrayForSparse);\r      else{\r	fHistMassPtImpParTCDs[2]->Fill(arrayForSparse);\r	fHistMassPtImpParTCDs[3]->Fill(arrayForSparseTrue);\r      }\r    }\r    \r    // fill PID\r    FillSpectrum(dstarD0pi,isDStar,fCuts,isSelected,fOutputPID);\r    SideBandBackground(dstarD0pi,fCuts,isSelected, fOutputPID);\r    //WrongSignForDStar(dstarD0pi,fCuts,fOutputPID);\r\r    //swich off the PID selection\r    fCuts->SetUsePID(kFALSE);\r    Int_t isSelectedNoPID=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected\r    fCuts->SetUsePID(kTRUE);\r\r    FillSpectrum(dstarD0pi,isDStar,fCuts,isSelectedNoPID,fOutputPID);\r    SideBandBackground(dstarD0pi,fCuts,isSelectedNoPID, fOutputPID);\r\r    // rare D search ------ \r    if(fDoSearch){\r      TLorentzVector LorentzTrack1(0,0,0,0); // lorentz 4 vector\r      TLorentzVector LorentzTrack2(0,0,0,0); // lorentz 4 vector\r      \r      for (Int_t i=0; i<aodEvent->GetNTracks(); i++){ \r	\r	AliAODTrack* aodTrack = aodEvent->GetTrack(i);\r	\r	if(dstarD0pi->Charge() == aodTrack->Charge()) continue;\r	if((!(aodTrack->GetStatus()&AliESDtrack::kITSrefit)|| (!(aodTrack->GetStatus()&AliESDtrack::kTPCrefit)))) continue;\r	if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))>0.02) continue;\r	\r	//build the D1 mass\r	Double_t mass = TDatabasePDG::Instance()->GetParticle(211)->Mass();\r\r	LorentzTrack1.SetPxPyPzE( dstarD0pi->Px(),dstarD0pi->Py(), dstarD0pi->Pz(), dstarD0pi->E(413) );\r	LorentzTrack2.SetPxPyPzE( aodTrack->Px(),aodTrack->Py(), aodTrack->Pz(),aodTrack->E(mass) );\r	\r	//D1 mass\r	Double_t d1mass = ((LorentzTrack1+LorentzTrack2).M());\r	//mass difference - at 0.4117 and 0.4566\r	fDeltaMassD1->Fill(d1mass-dstarD0pi->InvMassDstarKpipi());\r      }\r    }\r\r    if(isDStar == 1) { \r      fTrueDiff2->Fill(dstarD0pi->Pt(),dstarD0pi->DeltaInvMass());\r    }\r    \r  }\r  \r  fCounter->StoreCandidates(aodEvent,nSelectedProd,kTRUE);  \r  fCounter->StoreCandidates(aodEvent,nSelectedAna,kFALSE); \r\r  AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));\r  \r  PostData(1,fOutput);\r  PostData(2,fOutputAll);\r  PostData(3,fOutputPID);\r  PostData(5,fCounter);\r\r}\r//________________________________________ terminate ___________________________\rvoid AliAnalysisTaskSEDStarSpectra::Terminate(Option_t*)\r{    \r  // The Terminate() function is the last function to be called during\r  // a query. It always runs on the client, it can be used to present\r  // the results graphically or save the results to file.\r  \r  //Info("Terminate","");\r  AliAnalysisTaskSE::Terminate();\r  \r  fOutput = dynamic_cast<TList*> (GetOutputData(1));\r  if (!fOutput) {     \r    printf("ERROR: fOutput not available\n");\r    return;\r  }\r  \r  fCEvents        = dynamic_cast<TH1F*>(fOutput->FindObject("fCEvents"));\r  fDeltaMassD1     = dynamic_cast<TH1F*>(fOutput->FindObject("fDeltaMassD1"));\r  fTrueDiff2      = dynamic_cast<TH2F*>(fOutput->FindObject("fTrueDiff2"));\r\r  fOutputAll = dynamic_cast<TList*> (GetOutputData(1));\r  if (!fOutputAll) {\r    printf("ERROR: fOutputAll not available\n");\r    return;\r  }\r  fOutputPID = dynamic_cast<TList*> (GetOutputData(2));\r  if (!fOutputPID) {\r    printf("ERROR: fOutputPID not available\n");\r    return;\r  }\r \r\r  return;\r}\r//___________________________________________________________________________\rvoid AliAnalysisTaskSEDStarSpectra::UserCreateOutputObjects() { \r // output\r  Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());\r  \r  //slot #1  \r  //OpenFile(1);\r  fOutput = new TList();\r  fOutput->SetOwner();\r  fOutput->SetName("chist0");\r\r  fOutputAll = new TList();\r  fOutputAll->SetOwner();\r  fOutputAll->SetName("listAll");\r\r  fOutputPID = new TList();\r  fOutputPID->SetOwner();\r  fOutputPID->SetName("listPID");\r    \r  // define histograms\r  DefineHistograms();\r\r //Counter for Normalization\r fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));\r\r if(fDoImpParDstar) CreateImpactParameterHistos();\r\r  PostData(1,fOutput);\r  PostData(2,fOutputAll);\r  PostData(3,fOutputPID);\r\r  return;\r}\r//___________________________________ hiostograms _______________________________________\rvoid  AliAnalysisTaskSEDStarSpectra::DefineHistograms(){\r\r  fCEvents = new TH1F("fCEvents","conter",11,0,11);\r  fCEvents->SetStats(kTRUE);\r  fCEvents->GetXaxis()->SetTitle("1");\r  fCEvents->GetYaxis()->SetTitle("counts");\r  fOutput->Add(fCEvents);\r\r  fTrueDiff2 = new TH2F("DiffDstar_pt","True Reco diff vs pt",200,0,15,900,0,0.3);\r  fOutput->Add(fTrueDiff2);\r\r  fDeltaMassD1 = new TH1F("DeltaMassD1","delta mass d1",600,0,0.8);\r  fOutput->Add(fDeltaMassD1);\r\r  const Int_t nhist=14;\r  TString nameMass=" ", nameSgn=" ", nameBkg=" ";\r\r  for(Int_t i=-2;i<nhist;i++){\r    nameMass="histDeltaMass_";\r    nameMass+=i+1;\r    nameSgn="histDeltaSgn_";\r    nameSgn+=i+1;\r    nameBkg="histDeltaBkg_";\r    nameBkg+=i+1; \r    \r    if (i==-2) {\r      nameMass="histDeltaMass";\r      nameSgn="histDeltaSgn";\r      nameBkg="histDeltaBkg";\r    }\r    \r    TH1F* spectrumMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} invariant mass; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);\r    TH1F* spectrumSgn = new TH1F(nameSgn.Data(), "D^{*}-D^{0} Signal invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);\r    TH1F* spectrumBkg = new TH1F(nameBkg.Data(), "D^{*}-D^{0} Background invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);\r    \r    nameMass="histD0Mass_";\r    nameMass+=i+1;\r    nameSgn="histD0Sgn_";\r    nameSgn+=i+1;\r    nameBkg="histD0Bkg_";\r    nameBkg+=i+1;\r    \r    if (i==-2) {\r      nameMass="histD0Mass";\r      nameSgn="histD0Sgn";\r      nameBkg="histD0Bkg";\r    }\r\r    TH1F* spectrumD0Mass = new TH1F(nameMass.Data(),"D^{0} invariant mass; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);\r    TH1F* spectrumD0Sgn = new TH1F(nameSgn.Data(), "D^{0} Signal invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);\r    TH1F* spectrumD0Bkg = new TH1F(nameBkg.Data(), "D^{0} Background invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);\r\r    nameMass="histDstarMass_";\r    nameMass+=i+1;\r    nameSgn="histDstarSgn_";\r    nameSgn+=i+1;\r    nameBkg="histDstarBkg_";\r    nameBkg+=i+1;\r\r    if (i==-2) {\r      nameMass="histDstarMass";\r      nameSgn="histDstarSgn";\r      nameBkg="histDstarBkg";\r    }\r\r    TH1F* spectrumDstarMass = new TH1F(nameMass.Data(),"D^{*} invariant mass; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);\r    TH1F* spectrumDstarSgn = new TH1F(nameSgn.Data(), "D^{*} Signal invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);\r    TH1F* spectrumDstarBkg = new TH1F(nameBkg.Data(), "D^{*} Background invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);\r\r    nameMass="histSideBandMass_";\r    nameMass+=i+1;\r    if (i==-2) { \r      nameMass="histSideBandMass";\r    }\r    \r    TH1F* spectrumSideBandMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} sideband mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);\r\r    nameMass="histWrongSignMass_";\r    nameMass+=i+1;\r    if (i==-2) { \r      nameMass="histWrongSignMass";\r    }\r    \r    TH1F* spectrumWrongSignMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} wrongsign mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);\r\r\r    spectrumMass->Sumw2();\r    spectrumSgn->Sumw2();\r    spectrumBkg->Sumw2();\r    \r    spectrumMass->SetLineColor(6);\r    spectrumSgn->SetLineColor(2);\r    spectrumBkg->SetLineColor(4);\r    \r    spectrumMass->SetMarkerStyle(20);\r    spectrumSgn->SetMarkerStyle(20);\r    spectrumBkg->SetMarkerStyle(20);\r    spectrumMass->SetMarkerSize(0.6);\r    spectrumSgn->SetMarkerSize(0.6);\r    spectrumBkg->SetMarkerSize(0.6);\r    spectrumMass->SetMarkerColor(6);\r    spectrumSgn->SetMarkerColor(2);\r    spectrumBkg->SetMarkerColor(4);\r\r    spectrumD0Mass->Sumw2();\r    spectrumD0Sgn->Sumw2();\r    spectrumD0Bkg->Sumw2();\r\r    spectrumD0Mass->SetLineColor(6);\r    spectrumD0Sgn->SetLineColor(2);\r    spectrumD0Bkg->SetLineColor(4);\r\r    spectrumD0Mass->SetMarkerStyle(20);\r    spectrumD0Sgn->SetMarkerStyle(20);\r    spectrumD0Bkg->SetMarkerStyle(20);\r    spectrumD0Mass->SetMarkerSize(0.6);\r    spectrumD0Sgn->SetMarkerSize(0.6);\r    spectrumD0Bkg->SetMarkerSize(0.6);\r    spectrumD0Mass->SetMarkerColor(6);\r    spectrumD0Sgn->SetMarkerColor(2);\r    spectrumD0Bkg->SetMarkerColor(4);\r\r    spectrumDstarMass->Sumw2();\r    spectrumDstarSgn->Sumw2();\r    spectrumDstarBkg->Sumw2();\r\r    spectrumDstarMass->SetLineColor(6);\r    spectrumDstarSgn->SetLineColor(2);\r    spectrumDstarBkg->SetLineColor(4);\r\r    spectrumDstarMass->SetMarkerStyle(20);\r    spectrumDstarSgn->SetMarkerStyle(20);\r    spectrumDstarBkg->SetMarkerStyle(20);\r    spectrumDstarMass->SetMarkerSize(0.6);\r    spectrumDstarSgn->SetMarkerSize(0.6);\r    spectrumDstarBkg->SetMarkerSize(0.6);\r    spectrumDstarMass->SetMarkerColor(6);\r    spectrumDstarSgn->SetMarkerColor(2);\r    spectrumDstarBkg->SetMarkerColor(4);\r\r    spectrumSideBandMass->Sumw2();\r    spectrumSideBandMass->SetLineColor(4);\r    spectrumSideBandMass->SetMarkerStyle(20);\r    spectrumSideBandMass->SetMarkerSize(0.6);\r    spectrumSideBandMass->SetMarkerColor(4);\r\r    spectrumWrongSignMass->Sumw2();\r    spectrumWrongSignMass->SetLineColor(4);\r    spectrumWrongSignMass->SetMarkerStyle(20);\r    spectrumWrongSignMass->SetMarkerSize(0.6);\r    spectrumWrongSignMass->SetMarkerColor(4);\r\r    TH1F* allMass = (TH1F*)spectrumMass->Clone();\r    TH1F* allSgn  = (TH1F*)spectrumSgn->Clone();\r    TH1F* allBkg  = (TH1F*)spectrumBkg->Clone();\r\r    TH1F* pidMass = (TH1F*)spectrumMass->Clone();\r    TH1F* pidSgn  = (TH1F*)spectrumSgn->Clone();\r    TH1F* pidBkg  = (TH1F*)spectrumBkg->Clone();\r\r    fOutputAll->Add(allMass);\r    fOutputAll->Add(allSgn);\r    fOutputAll->Add(allBkg);\r\r    fOutputPID->Add(pidMass);\r    fOutputPID->Add(pidSgn);\r    fOutputPID->Add(pidBkg);\r\r    TH1F* allD0Mass = (TH1F*)spectrumD0Mass->Clone();\r    TH1F* allD0Sgn  = (TH1F*)spectrumD0Sgn->Clone();\r    TH1F* allD0Bkg  = (TH1F*)spectrumD0Bkg->Clone();\r\r    TH1F* pidD0Mass = (TH1F*)spectrumD0Mass->Clone();\r    TH1F* pidD0Sgn  = (TH1F*)spectrumD0Sgn->Clone();\r    TH1F* pidD0Bkg  = (TH1F*)spectrumD0Bkg->Clone();\r\r    fOutputAll->Add(allD0Mass);\r    fOutputAll->Add(allD0Sgn);\r    fOutputAll->Add(allD0Bkg);\r\r    fOutputPID->Add(pidD0Mass);\r    fOutputPID->Add(pidD0Sgn);\r    fOutputPID->Add(pidD0Bkg);\r  \r    TH1F* allDstarMass = (TH1F*)spectrumDstarMass->Clone();\r    TH1F* allDstarSgn = (TH1F*)spectrumDstarSgn->Clone();\r    TH1F* allDstarBkg = (TH1F*)spectrumDstarBkg->Clone();\r    \r    TH1F* pidDstarMass = (TH1F*)spectrumDstarMass->Clone();\r    TH1F* pidDstarSgn = (TH1F*)spectrumDstarSgn->Clone();\r    TH1F* pidDstarBkg = (TH1F*)spectrumDstarBkg->Clone();\r    \r    fOutputAll->Add(allDstarMass);\r    fOutputAll->Add(allDstarSgn);\r    fOutputAll->Add(allDstarBkg);\r\r    fOutputPID->Add(pidDstarMass);\r    fOutputPID->Add(pidDstarSgn);\r    fOutputPID->Add(pidDstarBkg);\r\r    TH1F* allSideBandMass = (TH1F*)spectrumSideBandMass->Clone();\r    TH1F* pidSideBandMass = (TH1F*)spectrumSideBandMass->Clone();\r\r    fOutputAll->Add(allSideBandMass);\r    fOutputPID->Add(pidSideBandMass);\r   \r    TH1F* allWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();\r    TH1F* pidWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();\r\r    fOutputAll->Add(allWrongSignMass);\r    fOutputPID->Add(pidWrongSignMass);\r   \r  }\r  \r  // pt spectra\r  nameMass="ptMass";\r  nameSgn="ptSgn";\r  nameBkg="ptBkg";\r  \r  TH1F* ptspectrumMass = new TH1F(nameMass.Data(),"D^{*} p_{T}; p_{T} [GeV]; Entries",200,0,10);\r  TH1F* ptspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal p_{T} - MC; p_{T} [GeV]; Entries",200,0,10);\r  TH1F* ptspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background p_{T} - MC; p_{T} [GeV]; Entries",200,0,10);\r  \r  ptspectrumMass->Sumw2();\r  ptspectrumSgn->Sumw2();\r  ptspectrumBkg->Sumw2();\r  \r  ptspectrumMass->SetLineColor(6);\r  ptspectrumSgn->SetLineColor(2);\r  ptspectrumBkg->SetLineColor(4);\r  \r  ptspectrumMass->SetMarkerStyle(20);\r  ptspectrumSgn->SetMarkerStyle(20);\r  ptspectrumBkg->SetMarkerStyle(20);\r  ptspectrumMass->SetMarkerSize(0.6);\r  ptspectrumSgn->SetMarkerSize(0.6);\r  ptspectrumBkg->SetMarkerSize(0.6);\r  ptspectrumMass->SetMarkerColor(6);\r  ptspectrumSgn->SetMarkerColor(2);\r  ptspectrumBkg->SetMarkerColor(4);\r  \r  TH1F* ptallMass = (TH1F*)ptspectrumMass->Clone();\r  TH1F* ptallSgn = (TH1F*)ptspectrumSgn->Clone();\r  TH1F* ptallBkg = (TH1F*)ptspectrumBkg->Clone();\r  \r  TH1F* ptpidMass = (TH1F*)ptspectrumMass->Clone();\r  TH1F* ptpidSgn = (TH1F*)ptspectrumSgn->Clone();\r  TH1F* ptpidBkg = (TH1F*)ptspectrumBkg->Clone();\r    \r  fOutputAll->Add(ptallMass);\r  fOutputAll->Add(ptallSgn);\r  fOutputAll->Add(ptallBkg);\r  \r  fOutputPID->Add(ptpidMass);\r  fOutputPID->Add(ptpidSgn);\r  fOutputPID->Add(ptpidBkg);\r  \r  // eta spectra\r  nameMass="etaMass";\r  nameSgn="etaSgn";\r  nameBkg="etaBkg";\r  \r  TH1F* etaspectrumMass = new TH1F(nameMass.Data(),"D^{*} #eta; #eta; Entries",200,-1,1);\r  TH1F* etaspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal #eta - MC; #eta; Entries",200,-1,1);\r  TH1F* etaspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background #eta - MC; #eta; Entries",200,-1,1);\r  \r  etaspectrumMass->Sumw2();\r  etaspectrumSgn->Sumw2();\r  etaspectrumBkg->Sumw2();\r  \r  etaspectrumMass->SetLineColor(6);\r  etaspectrumSgn->SetLineColor(2);\r  etaspectrumBkg->SetLineColor(4);\r  \r  etaspectrumMass->SetMarkerStyle(20);\r  etaspectrumSgn->SetMarkerStyle(20);\r  etaspectrumBkg->SetMarkerStyle(20);\r  etaspectrumMass->SetMarkerSize(0.6);\r  etaspectrumSgn->SetMarkerSize(0.6);\r  etaspectrumBkg->SetMarkerSize(0.6);\r  etaspectrumMass->SetMarkerColor(6);\r  etaspectrumSgn->SetMarkerColor(2);\r  etaspectrumBkg->SetMarkerColor(4);\r  \r  TH1F* etaallMass = (TH1F*)etaspectrumMass->Clone();\r  TH1F* etaallSgn = (TH1F*)etaspectrumSgn->Clone();\r  TH1F* etaallBkg = (TH1F*)etaspectrumBkg->Clone();\r  \r  TH1F* etapidMass = (TH1F*)etaspectrumMass->Clone();\r  TH1F* etapidSgn = (TH1F*)etaspectrumSgn->Clone();\r  TH1F* etapidBkg = (TH1F*)etaspectrumBkg->Clone();\r    \r  fOutputAll->Add(etaallMass);\r  fOutputAll->Add(etaallSgn);\r  fOutputAll->Add(etaallBkg);\r  \r  fOutputPID->Add(etapidMass);\r  fOutputPID->Add(etapidSgn);\r  fOutputPID->Add(etapidBkg);\r  \r  return;\r}\r//________________________________________________________________________\rvoid AliAnalysisTaskSEDStarSpectra::FillSpectrum(AliAODRecoCascadeHF *part, Int_t isDStar, AliRDHFCutsDStartoKpipi *cuts,Int_t isSel, TList *listout){\r  //\r  // Fill histos for D* spectrum\r  //\r  \r  if(!isSel) return;\r\r  // D0 window\r  Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();\r  Double_t invmassD0   = part->InvMassD0();  \r  if (TMath::Abs(invmassD0-mPDGD0)>fD0Window) return; \r\r\r  Int_t ptbin=cuts->PtBin(part->Pt());  \r  Double_t pt = part->Pt();\r  Double_t eta = part->Eta();  \r  \r  Double_t invmassDelta = part->DeltaInvMass();\r  Double_t invmassDstar = part->InvMassDstarKpipi();\r  \r  TString fillthis="";\r  Bool_t massInRange=kFALSE;\r  \r  Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();\r  \r  // delta M(Kpipi)-M(Kpi)\r  if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<fPeakWindow) massInRange=kTRUE;  \r  \r  if(fUseMCInfo) {\r    if(isDStar==1) {\r      fillthis="histD0Sgn_";\r      fillthis+=ptbin;\r      ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);\r      fillthis="histD0Sgn";\r      ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);\r      fillthis="histDstarSgn_";\r      fillthis+=ptbin;\r      ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);\r      fillthis="histDstarSgn";\r      ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);\r      fillthis="histDeltaSgn_";\r      fillthis+=ptbin;\r      ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);\r      fillthis="histDeltaSgn";\r      ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);\r    if (massInRange) {\r	fillthis="ptSgn";\r	((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);\r	fillthis="etaSgn";\r	((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);\r      }\r    }\r    else {//background\r      fillthis="histD0Bkg_";\r      fillthis+=ptbin;\r      ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);\r      fillthis="histD0Bkg";\r      ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);\r      fillthis="histDstarBkg_";\r      fillthis+=ptbin;\r      ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);\r      fillthis="histDstarBkg";\r      ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);\r      fillthis="histDeltaBkg_";\r      fillthis+=ptbin;\r      ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);\r      fillthis="histDeltaBkg";\r      ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);\r     if (massInRange) {\r	fillthis="ptBkg";\r	((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);\r	fillthis="etaBkg";\r	((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);\r      }\r    }\r  }\r  //no MC info, just cut selection\r  fillthis="histD0Mass_";\r  fillthis+=ptbin;\r  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);\r  fillthis="histD0Mass";\r  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);\r  fillthis="histDstarMass_";\r  fillthis+=ptbin;\r  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);\r  fillthis="histDstarMass";\r  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);\r  fillthis="histDeltaMass_";\r  fillthis+=ptbin;\r  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);\r  fillthis="histDeltaMass";\r  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);\r  \r  if (massInRange) {\r    fillthis="ptMass";\r    ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);\r    fillthis="etaMass";\r    ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);\r  }\r \r  return;\r}\r//______________________________ side band background for D*___________________________________\rvoid AliAnalysisTaskSEDStarSpectra::SideBandBackground(AliAODRecoCascadeHF *part,  AliRDHFCutsDStartoKpipi *cuts, Int_t isSel, TList *listout){\r\r  //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas \r  // (expected detector resolution) on the left and right frm the D0 mass. Each band\r  //  has a width of ~5 sigmas. Two band needed  for opening angle considerations   \r\r  if(!isSel) return;\r\r  Int_t ptbin=cuts->PtBin(part->Pt());\r  \r  Bool_t massInRange=kFALSE;\r  \r  // select the side bands intervall\r  Double_t invmassD0    = part->InvMassD0();\r  if(TMath::Abs(invmassD0-1.865)>4*fD0Window && TMath::Abs(invmassD0-1.865)<8*fD0Window){\r    \r    // for pt and eta\r    Double_t invmassDelta = part->DeltaInvMass();\r    if (TMath::Abs(invmassDelta-0.14557)<fPeakWindow) massInRange=kTRUE;\r    \r    TString fillthis="";\r    fillthis="histSideBandMass_";\r    fillthis+=ptbin;\r    ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);\r    fillthis="histSideBandMass";\r    ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);\r    \r  }\r}\r//________________________________________________________________________________________________________________\rvoid AliAnalysisTaskSEDStarSpectra::WrongSignForDStar(AliAODRecoCascadeHF *part,  AliRDHFCutsDStartoKpipi *cuts, TList *listout){\r  //\r  // assign the wrong charge to the soft pion to create background\r  //\r  Int_t ptbin=cuts->PtBin(part->Pt());\r  \r  Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();\r  Double_t invmassD0   = part->InvMassD0();  \r  if (TMath::Abs(invmassD0-mPDGD0)>fD0Window) return; \r\r  AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)part->Get2Prong();\r\r  Int_t okD0WrongSign,okD0barWrongSign;\r  Double_t wrongMassD0=0.;\r  \r  Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kCandidate); //selected\r   if (!isSelected){\r    return;\r  }\r\r  okD0WrongSign =  1;\r  okD0barWrongSign = 1;\r  \r  //if is D*+ than assume D0bar\r  if(part->Charge()>0 && (isSelected ==1)) { \r    okD0WrongSign = 0;\r  }\r  if(part->Charge()<0 && (isSelected ==2)){\r    okD0barWrongSign = 0;\r  }\r  \r  // assign the wrong mass in case the cuts return both D0 and D0bar\r  if(part->Charge()>0 && (isSelected ==3)) { \r    okD0WrongSign = 0;\r  } else if(part->Charge()<0 && (isSelected ==3)){\r    okD0barWrongSign = 0;\r  }\r  \r  //wrong D0 inv mass\r  if(okD0WrongSign!=0){\r    wrongMassD0 = theD0particle->InvMassD0();\r  }else if(okD0WrongSign==0){\r    wrongMassD0 = theD0particle->InvMassD0bar();\r  }\r  \r  if(TMath::Abs(wrongMassD0-1.865)<fD0Window){\r    \r    // wrong D* inv mass   \r    Double_t e[3];\r    if (part->Charge()>0){\r      e[0]=theD0particle->EProng(0,321);\r      e[1]=theD0particle->EProng(1,211);\r    }else{\r      e[0]=theD0particle->EProng(0,211);\r      e[1]=theD0particle->EProng(1,321);\r    }\r    e[2]=part->EProng(0,211);\r    \r    Double_t esum = e[0]+e[1]+e[2];\r    Double_t pds = part->P();\r\r    Double_t   wrongMassDstar = TMath::Sqrt(esum*esum-pds*pds);\r\r    TString fillthis="";\r    fillthis="histWrongSignMass_";\r    fillthis+=ptbin;\r    ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);\r    fillthis="histWrongSignMass";\r    ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);\r    \r  }\r}\r//-------------------------------------------------------------------------------\rInt_t AliAnalysisTaskSEDStarSpectra::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {		\r  //\r  // checking whether the mother of the particles come from a charm or a bottom quark\r  //\r	\r  Int_t pdgGranma = 0;\r  Int_t mother = 0;\r  mother = mcPartCandidate->GetMother();\r  Int_t istep = 0;\r  Int_t abspdgGranma =0;\r  Bool_t isFromB=kFALSE;\r  Bool_t isQuarkFound=kFALSE;\r  while (mother >0 ){\r    istep++;\r    AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));\r    if (mcGranma){\r      pdgGranma = mcGranma->GetPdgCode();\r      abspdgGranma = TMath::Abs(pdgGranma);\r      if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){\r	isFromB=kTRUE;\r      }\r      if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;\r      mother = mcGranma->GetMother();\r    }else{\r      AliError("Failed casting the mother particle!");\r      break;\r    }\r  }\r  \r  if(isFromB) return 5;\r  else return 4;\r}\r//-------------------------------------------------------------------------------------\rFloat_t AliAnalysisTaskSEDStarSpectra::GetTrueImpactParameterD0(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partDp) const {\r  // true impact parameter calculation\r\r  Double_t vtxTrue[3];\r  mcHeader->GetVertex(vtxTrue);\r  Double_t origD[3];\r  partDp->XvYvZv(origD);	  \r  Short_t charge=partDp->Charge();\r  Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];\r  Int_t labelFirstDau = partDp->GetDaughter(0); \r\r  Int_t nDau=partDp->GetNDaughters();\r\r  Int_t theDau=0;\r  if(nDau==2){\r    for(Int_t iDau=0; iDau<2; iDau++){\r      Int_t ind = labelFirstDau+iDau;\r      AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));\r      if(!part){\r	AliError("Daughter particle not found in MC array");\r	return 99999.;\r      } \r      Int_t pdgCode=TMath::Abs(part->GetPdgCode());\r      if(pdgCode==211 || pdgCode==321){\r	pXdauTrue[theDau]=part->Px();\r	pYdauTrue[theDau]=part->Py();\r	pZdauTrue[theDau]=part->Pz();\r	++theDau;\r      }\r    }\r  }\r  if(theDau!=2){\r    AliError("Wrong number of decay prongs");\r    return 99999.;\r  }\r\r  Double_t d0dummy[3]={0.,0.,0.};\r  AliAODRecoDecayHF aodD0MC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);\r  return aodD0MC.ImpParXY();\r\r}\r//______________________________________________________-\rvoid AliAnalysisTaskSEDStarSpectra::CreateImpactParameterHistos(){\r  // Histos for impact paramter study\r\r  Int_t nbins[3]={400,200,fNImpParBins};\r  Double_t xmin[3]={1.75,0.,fLowerImpPar};\r  Double_t xmax[3]={1.98,20.,fHigherImpPar};\r\r  fHistMassPtImpParTCDs[0]=new THnSparseF("hMassPtImpParAll",\r					"Mass vs. pt vs.imppar - All",\r					3,nbins,xmin,xmax);\r  fHistMassPtImpParTCDs[1]=new THnSparseF("hMassPtImpParPrompt",\r					"Mass vs. pt vs.imppar - promptD",\r					3,nbins,xmin,xmax);\r  fHistMassPtImpParTCDs[2]=new THnSparseF("hMassPtImpParBfeed",\r					"Mass vs. pt vs.imppar - DfromB",\r					3,nbins,xmin,xmax);\r  fHistMassPtImpParTCDs[3]=new THnSparseF("hMassPtImpParTrueBfeed",\r					"Mass vs. pt vs.true imppar -DfromB",\r					3,nbins,xmin,xmax);\r  fHistMassPtImpParTCDs[4]=new THnSparseF("hMassPtImpParBkg",\r				        "Mass vs. pt vs.imppar - backgr.",\r					3,nbins,xmin,xmax);\r\r  for(Int_t i=0; i<5;i++){\r    fOutput->Add(fHistMassPtImpParTCDs[i]);\r  }\r}\r