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