-/**************************************************************************\r
- * Copyright(c) 2007-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
- * appear 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
-// Analysis task to produce Ds candidates mass spectra //\r
-// Origin: F.Prino, Torino, prino@to.infn.it //\r
-// //\r
-///////////////////////////////////////////////////////////////////\r
-\r
-#include <TClonesArray.h>\r
-#include <TNtuple.h>\r
-#include <TList.h>\r
-#include <TString.h>\r
-#include <TH1F.h>\r
-#include <TMath.h>\r
-#include <TDatabasePDG.h>\r
-\r
-#include "AliAnalysisManager.h"\r
-#include "AliAODHandler.h"\r
-#include "AliAODEvent.h"\r
-#include "AliAODVertex.h"\r
-#include "AliAODTrack.h"\r
-#include "AliAODMCHeader.h"\r
-#include "AliAODMCParticle.h"\r
-#include "AliAODRecoDecayHF3Prong.h"\r
-#include "AliAnalysisVertexingHF.h"\r
-#include "AliRDHFCutsDstoKKpi.h"\r
-#include "AliAnalysisTaskSE.h"\r
-#include "AliNormalizationCounter.h"\r
-#include "AliAnalysisTaskSEDs.h"\r
-\r
-ClassImp(AliAnalysisTaskSEDs)\r
-\r
-\r
-//________________________________________________________________________\r
-AliAnalysisTaskSEDs::AliAnalysisTaskSEDs():\r
- AliAnalysisTaskSE(),\r
- fOutput(0), \r
- fHistNEvents(0),\r
- fPtVsMass(0),\r
- fYVsPt(0),\r
- fYVsPtSig(0),\r
- fNtupleDs(0),\r
- fFillNtuple(0),\r
- fReadMC(kFALSE),\r
- fDoCutVarHistos(kTRUE),\r
- fUseSelectionBit(kFALSE),\r
- fNPtBins(0),\r
- fListCuts(0),\r
- fMassRange(0.8),\r
- fMassBinSize(0.002),\r
- fCounter(0),\r
- fProdCuts(0),\r
- fAnalysisCuts(0)\r
-{\r
- // Default constructor\r
-\r
- for(Int_t i=0;i<4;i++) {\r
- if(fChanHist[i]) fChanHist[i]=0;\r
- }\r
-\r
- for(Int_t i=0;i<4*kMaxPtBins;i++){\r
- \r
- if(fPtCandHist[i]) fPtCandHist[i]=0;\r
- if(fMassHist[i]) fMassHist[i]=0;\r
- if(fCosPHist[i]) fCosPHist[i]=0;\r
- if(fDLenHist[i]) fDLenHist[i]=0;\r
- if(fSumd02Hist[i]) fSumd02Hist[i]=0;\r
- if(fSigVertHist[i]) fSigVertHist[i]=0;\r
- if(fPtMaxHist[i]) fPtMaxHist[i]=0;\r
- if(fDCAHist[i]) fDCAHist[i]=0;\r
- if(fPtProng0Hist[i]) fPtProng0Hist[i]=0;\r
- if(fPtProng1Hist[i]) fPtProng1Hist[i]=0;\r
- if(fPtProng2Hist[i]) fPtProng2Hist[i]=0;\r
- if(fDalitz[i]) fDalitz[i]=0;\r
- if(fDalitzPhi[i]) fDalitzPhi[i]=0;\r
- if(fDalitzK0st[i]) fDalitzK0st[i]=0;\r
-\r
- }\r
- for(Int_t i=0;i<3*kMaxPtBins;i++){\r
- if(fMassHistPhi[i]) fMassHistPhi[i]=0;\r
- if(fMassHistK0st[i]) fMassHistK0st[i]=0;\r
- \r
- }\r
-\r
- for(Int_t i=0;i<kMaxPtBins+1;i++){\r
- fPtLimits[i]=0;\r
- }\r
-\r
-}\r
-\r
-//________________________________________________________________________\r
-AliAnalysisTaskSEDs::AliAnalysisTaskSEDs(const char *name, AliRDHFCutsDstoKKpi* productioncuts, AliRDHFCutsDstoKKpi* analysiscuts,Int_t fillNtuple):\r
- AliAnalysisTaskSE(name),\r
- fOutput(0),\r
- fHistNEvents(0),\r
- fPtVsMass(0),\r
- fYVsPt(0),\r
- fYVsPtSig(0),\r
- fNtupleDs(0),\r
- fFillNtuple(fillNtuple),\r
- fReadMC(kFALSE),\r
- fDoCutVarHistos(kTRUE),\r
- fUseSelectionBit(kFALSE),\r
- fNPtBins(0),\r
- fListCuts(0),\r
- fMassRange(0.8),\r
- fMassBinSize(0.002),\r
- fCounter(0),\r
- fProdCuts(productioncuts),\r
- fAnalysisCuts(analysiscuts)\r
-{\r
- // Default constructor\r
- // Output slot #1 writes into a TList container\r
- \r
- for(Int_t i=0;i<4;i++) {\r
- if(fChanHist[i]) fChanHist[i]=0;\r
- }\r
-\r
- for(Int_t i=0;i<4*kMaxPtBins;i++){\r
- \r
- if(fPtCandHist[i]) fPtCandHist[i]=0;\r
- if(fMassHist[i]) fMassHist[i]=0;\r
- if(fCosPHist[i]) fCosPHist[i]=0;\r
- if(fDLenHist[i]) fDLenHist[i]=0;\r
- if(fSumd02Hist[i]) fSumd02Hist[i]=0;\r
- if(fSigVertHist[i]) fSigVertHist[i]=0;\r
- if(fPtMaxHist[i]) fPtMaxHist[i]=0;\r
- if(fDCAHist[i]) fDCAHist[i]=0;\r
- if(fPtProng0Hist[i]) fPtProng0Hist[i]=0;\r
- if(fPtProng1Hist[i]) fPtProng1Hist[i]=0;\r
- if(fPtProng2Hist[i]) fPtProng2Hist[i]=0;\r
- if(fDalitz[i]) fDalitz[i]=0;\r
- if(fDalitzPhi[i]) fDalitzPhi[i]=0;\r
- if(fDalitzK0st[i]) fDalitzK0st[i]=0;\r
-\r
- }\r
- for(Int_t i=0;i<3*kMaxPtBins;i++){\r
- if(fMassHistPhi[i]) fMassHistPhi[i]=0;\r
- if(fMassHistK0st[i]) fMassHistK0st[i]=0;\r
- \r
- }\r
-\r
- for(Int_t i=0;i<kMaxPtBins+1;i++){\r
- fPtLimits[i]=0;\r
- }\r
- \r
- Int_t nptbins=fAnalysisCuts->GetNPtBins();\r
- Float_t *ptlim=fAnalysisCuts->GetPtBinLimits();\r
- SetPtBins(nptbins,ptlim);\r
-\r
- DefineOutput(1,TList::Class()); //My private output\r
-\r
- DefineOutput(2,TList::Class());\r
-\r
- DefineOutput(3,AliNormalizationCounter::Class());\r
- \r
- if(fFillNtuple>0){\r
- // Output slot #4 writes into a TNtuple container\r
- DefineOutput(4,TNtuple::Class()); //My private output\r
- }\r
- \r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskSEDs::SetPtBins(Int_t n, Float_t* lim){\r
- // define pt bins for analysis\r
- if(n>kMaxPtBins){\r
- printf("Max. number of Pt bins = %d\n",kMaxPtBins);\r
- fNPtBins=kMaxPtBins;\r
- fPtLimits[0]=0.;\r
- fPtLimits[1]=1.;\r
- fPtLimits[2]=3.;\r
- fPtLimits[3]=5.;\r
- fPtLimits[4]=10.;\r
- for(Int_t i=5; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;\r
- }else{\r
- fNPtBins=n;\r
- for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];\r
- for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;\r
- }\r
- if(fDebug > 1){\r
- printf("Number of Pt bins = %d\n",fNPtBins);\r
- for(Int_t i=0; i<fNPtBins; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fPtLimits[i],fPtLimits[i+1]); \r
- }\r
-}\r
-//________________________________________________________________________\r
-AliAnalysisTaskSEDs::~AliAnalysisTaskSEDs()\r
-{\r
- // Destructor\r
- delete fOutput;\r
- delete fHistNEvents;\r
- delete fListCuts;\r
-\r
- for(Int_t i=0;i<4*fNPtBins;i++){\r
- \r
- if(fMassHist[i]){ delete fMassHist[i]; fMassHist[i]=0;}\r
- if(fMassHistPhi[i]){ delete fMassHistPhi[i]; fMassHistPhi[i]=0;}\r
- if(fMassHistK0st[i]){ delete fMassHistK0st[i]; fMassHistK0st[i]=0;}\r
- if(fCosPHist[i]){ delete fCosPHist[i]; fCosPHist[i]=0;}\r
- if(fDLenHist[i]){ delete fDLenHist[i]; fDLenHist[i]=0;}\r
- if(fSumd02Hist[i]){ delete fSumd02Hist[i]; fSumd02Hist[i]=0;}\r
- if(fSigVertHist[i]){ delete fSigVertHist[i]; fSigVertHist[i]=0;}\r
- if(fPtMaxHist[i]){ delete fPtMaxHist[i]; fPtMaxHist[i]=0;}\r
- if(fDCAHist[i]){ delete fDCAHist[i]; fDCAHist[i]=0;}\r
- if(fPtProng0Hist[i]){ delete fPtProng0Hist[i]; fPtProng0Hist[i]=0;}\r
- if(fPtProng1Hist[i]){ delete fPtProng1Hist[i]; fPtProng1Hist[i]=0;}\r
- if(fPtProng2Hist[i]){ delete fPtProng2Hist[i]; fPtProng2Hist[i]=0;}\r
- if(fDalitz[i]){ delete fDalitz[i]; fDalitz[i]=0;}\r
- if(fDalitzPhi[i]){ delete fDalitzPhi[i]; fDalitzPhi[i]=0;}\r
- if(fDalitzK0st[i]){ delete fDalitzK0st[i]; fDalitzK0st[i]=0;}\r
-\r
- }\r
-\r
- delete fPtVsMass;\r
- delete fYVsPt;\r
- delete fYVsPtSig;\r
- delete fNtupleDs;\r
- delete fCounter;\r
- delete fProdCuts;\r
- delete fAnalysisCuts;\r
-\r
-} \r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskSEDs::Init()\r
-{\r
- // Initialization\r
-\r
- if(fDebug > 1) printf("AnalysisTaskSEDs::Init() \n");\r
-\r
- fListCuts=new TList();\r
- fListCuts->SetOwner();\r
- fListCuts->SetName("CutObjects");\r
-\r
- AliRDHFCutsDstoKKpi *production = new AliRDHFCutsDstoKKpi(*fProdCuts);\r
- production->SetName("ProductionCuts");\r
- AliRDHFCutsDstoKKpi *analysis = new AliRDHFCutsDstoKKpi(*fAnalysisCuts);\r
- analysis->SetName("AnalysisCuts");\r
- \r
- fListCuts->Add(production);\r
- fListCuts->Add(analysis);\r
- PostData(2,fListCuts);\r
- return;\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskSEDs::UserCreateOutputObjects()\r
-{\r
- // Create the output container\r
- //\r
- if(fDebug > 1) printf("AnalysisTaskSEDs::UserCreateOutputObjects() \n");\r
-\r
- // Several histograms are more conveniently managed in a TList\r
- fOutput = new TList();\r
- fOutput->SetOwner();\r
- fOutput->SetName("OutputHistos");\r
-\r
- Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();\r
- \r
- Int_t nInvMassBins=(Int_t)(fMassRange/fMassBinSize+0.5);\r
- if(nInvMassBins%2==1) nInvMassBins++;\r
- // Double_t minMass=massDs-1.0*nInvMassBins*fMassBinSize;\r
- Double_t minMass=massDs-0.5*nInvMassBins*fMassBinSize;\r
- // Double_t maxMass=massDs+1.0*nInvMassBins*fMassBinSize;\r
- Double_t maxMass=massDs+0.5*nInvMassBins*fMassBinSize;\r
-\r
- TString hisname;\r
- TString htype;\r
- Int_t index;\r
- for(Int_t iType=0; iType<4; iType++){\r
- for(Int_t i=0;i<fNPtBins;i++){\r
- if(iType==0){\r
- htype="All";\r
- index=GetHistoIndex(i);\r
- }else if(iType==1){ \r
- htype="Sig";\r
- index=GetSignalHistoIndex(i);\r
- }else if(iType==2){ \r
- htype="Bkg";\r
- index=GetBackgroundHistoIndex(i);\r
- }else{ \r
- htype="ReflSig";\r
- index=GetReflSignalHistoIndex(i);\r
- }\r
- hisname.Form("hMass%sPt%d",htype.Data(),i);\r
- fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);\r
- fMassHist[index]->Sumw2();\r
- hisname.Form("hMass%sPt%dphi",htype.Data(),i);\r
- fMassHistPhi[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);\r
- fMassHistPhi[index]->Sumw2();\r
- hisname.Form("hMass%sPt%dk0st",htype.Data(),i);\r
- fMassHistK0st[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);\r
- fMassHistK0st[index]->Sumw2();\r
- hisname.Form("hCosP%sPt%d",htype.Data(),i);\r
- fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);\r
- fCosPHist[index]->Sumw2();\r
- hisname.Form("hDLen%sPt%d",htype.Data(),i);\r
- fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);\r
- fDLenHist[index]->Sumw2();\r
- hisname.Form("hSumd02%sPt%d",htype.Data(),i);\r
- fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);\r
- fSumd02Hist[index]->Sumw2();\r
- hisname.Form("hSigVert%sPt%d",htype.Data(),i);\r
- fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);\r
- fSigVertHist[index]->Sumw2();\r
- hisname.Form("hPtMax%sPt%d",htype.Data(),i);\r
- fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);\r
- fPtMaxHist[index]->Sumw2();\r
- hisname.Form("hPtCand%sPt%d",htype.Data(),i);\r
- fPtCandHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);\r
- fPtCandHist[index]->Sumw2();\r
- hisname.Form("hDCA%sPt%d",htype.Data(),i);\r
- fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);\r
- fDCAHist[index]->Sumw2();\r
- hisname.Form("hPtProng0%sPt%d",htype.Data(),i);\r
- fPtProng0Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);\r
- fPtProng0Hist[index]->Sumw2();\r
- hisname.Form("hPtProng1%sPt%d",htype.Data(),i);\r
- fPtProng1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);\r
- fPtProng1Hist[index]->Sumw2();\r
- hisname.Form("hPtProng2%sPt%d",htype.Data(),i);\r
- fPtProng2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);\r
- fPtProng2Hist[index]->Sumw2();\r
- hisname.Form("hDalitz%sPt%d",htype.Data(),i);\r
- fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);\r
- fDalitz[index]->Sumw2();\r
- hisname.Form("hDalitz%sPt%dphi",htype.Data(),i);\r
- fDalitzPhi[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);\r
- fDalitzPhi[index]->Sumw2();\r
- hisname.Form("hDalitz%sPt%dk0st",htype.Data(),i);\r
- fDalitzK0st[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);\r
- fDalitzK0st[index]->Sumw2();\r
- }\r
- }\r
-\r
- for(Int_t i=0; i<4*fNPtBins; i++){\r
- fOutput->Add(fMassHist[i]);\r
- fOutput->Add(fMassHistPhi[i]);\r
- fOutput->Add(fMassHistK0st[i]);\r
- fOutput->Add(fCosPHist[i]);\r
- fOutput->Add(fDLenHist[i]);\r
- fOutput->Add(fSumd02Hist[i]);\r
- fOutput->Add(fSigVertHist[i]);\r
- fOutput->Add(fPtMaxHist[i]);\r
- fOutput->Add(fPtCandHist[i]);\r
- fOutput->Add(fDCAHist[i]);\r
- fOutput->Add(fPtProng0Hist[i]);\r
- fOutput->Add(fPtProng1Hist[i]);\r
- fOutput->Add(fPtProng2Hist[i]);\r
- fOutput->Add(fDalitz[i]);\r
- fOutput->Add(fDalitzPhi[i]);\r
- fOutput->Add(fDalitzK0st[i]);\r
- }\r
-\r
- fChanHist[0] = new TH1F("hChanAll", "KKpi and piKK candidates",64,-0.5,63.5);\r
- fChanHist[1] = new TH1F("hChanSig", "KKpi and piKK candidates",64,-0.5,63.5);\r
- fChanHist[2] = new TH1F("hChanBkg", "KKpi and piKK candidates",64,-0.5,63.5);\r
- fChanHist[3] = new TH1F("hChanReflSig", "KKpi and piKK candidates",64,-0.5,63.5);\r
- for(Int_t i=0;i<4;i++){\r
- fChanHist[i]->Sumw2();\r
- fChanHist[i]->SetMinimum(0);\r
- fOutput->Add(fChanHist[i]);\r
- }\r
-\r
- fHistNEvents = new TH1F("hNEvents", "number of events ",12,-0.5,11.5);\r
- fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");\r
- fHistNEvents->GetXaxis()->SetBinLabel(2,"n. passing IsEvSelected");\r
- fHistNEvents->GetXaxis()->SetBinLabel(3,"n. rejected due to trigger");\r
- fHistNEvents->GetXaxis()->SetBinLabel(4,"n. rejected due to not reco vertex");\r
- fHistNEvents->GetXaxis()->SetBinLabel(5,"n. rejected for contr vertex");\r
- fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected for vertex out of accept");\r
- fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected for pileup events");\r
- fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events");\r
- fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of candidate");\r
- fHistNEvents->GetXaxis()->SetBinLabel(10,"no. of Ds after loose cuts");\r
- fHistNEvents->GetXaxis()->SetBinLabel(11,"no. of Ds after tight cuts");\r
- fHistNEvents->GetXaxis()->SetBinLabel(12,"no. of cand wo bitmask");\r
-\r
- fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);\r
-\r
- fHistNEvents->Sumw2();\r
- fHistNEvents->SetMinimum(0);\r
- fOutput->Add(fHistNEvents);\r
-\r
- fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nInvMassBins,minMass,maxMass,40,0.,20.);\r
- fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);\r
- fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);\r
-\r
- fOutput->Add(fPtVsMass);\r
- fOutput->Add(fYVsPt);\r
- fOutput->Add(fYVsPtSig);\r
-\r
- //Counter for Normalization\r
- fCounter = new AliNormalizationCounter("NormalizationCounter");\r
- fCounter->Init();\r
-\r
- PostData(1,fOutput); \r
- PostData(3,fCounter); \r
- \r
- if(fFillNtuple>0){\r
- OpenFile(4); // 4 is the slot number of the ntuple\r
- \r
- fNtupleDs = new TNtuple("fNtupleDs","Ds","labDs:retcode:pdgcode0:Pt0:Pt1:Pt2:PtRec:P0:P1:P2:PidTrackBit0:PidTrackBit1:PidTrackBit2:PointingAngle:PointingAngleXY:DecLeng:DecLengXY:NorDecLeng:NorDecLengXY:InvMassKKpi:InvMasspiKK:sigvert:d00:d01:d02:dca:d0square:InvMassPhiKKpi:InvMassPhipiKK:InvMassK0starKKpi:InvMassK0starpiKK:cosinePiDsFrameKKpi:cosinePiDsFramepiKK:cosineKPhiFrameKKpi:cosineKPhiFramepiKK:centrality"); \r
- \r
- }\r
- \r
-\r
- return;\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskSEDs::UserExec(Option_t */*option*/)\r
-{\r
- // Ds selection for current event, fill mass histos and selecetion variable histo\r
- // separate signal and backgound if fReadMC is activated\r
-\r
- AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());\r
-\r
- TClonesArray *array3Prong = 0;\r
- if(!aod && 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
- aod = 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
- array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");\r
- }\r
- } else if(aod) {\r
- array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");\r
- }\r
-\r
- if(!aod || !array3Prong) {\r
- printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");\r
- return;\r
- }\r
-\r
-\r
- // fix for temporary bug in ESDfilter \r
- // the AODs with null vertex pointer didn't pass the PhysSel\r
- if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;\r
-\r
-\r
- fHistNEvents->Fill(0); // count event\r
- // Post the data already here\r
- PostData(1,fOutput);\r
- \r
- fCounter->StoreEvent(aod,fProdCuts,fReadMC);\r
- \r
-\r
- Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);\r
- if(fAnalysisCuts->IsEventRejectedDueToTrigger())fHistNEvents->Fill(2);\r
- if(fAnalysisCuts->IsEventRejectedDueToNotRecoVertex())fHistNEvents->Fill(3);\r
- if(fAnalysisCuts->IsEventRejectedDueToVertexContributors())fHistNEvents->Fill(4);\r
- if(fAnalysisCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion())fHistNEvents->Fill(5);\r
- if(fAnalysisCuts->IsEventRejectedDueToPileupSPD())fHistNEvents->Fill(6);\r
- if(fAnalysisCuts->IsEventRejectedDueToCentrality())fHistNEvents->Fill(7);\r
- \r
- Float_t centrality=fAnalysisCuts->GetCentrality(aod);\r
- //Int_t runNumber=aod->GetRunNumber();\r
-\r
- if(!isEvSel)return;\r
- \r
- \r
-\r
- \r
- fHistNEvents->Fill(1);\r
-\r
- TClonesArray *arrayMC=0;\r
- AliAODMCHeader *mcHeader=0;\r
-\r
- // AOD primary vertex\r
- AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();\r
- // vtx1->Print();\r
- \r
- // load MC particles\r
- if(fReadMC){\r
- \r
- arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
- if(!arrayMC) {\r
- printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");\r
- return;\r
- }\r
- \r
- // load MC header\r
- mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());\r
- if(!mcHeader) {\r
- printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");\r
- return;\r
- }\r
- }\r
- \r
- Int_t n3Prong = array3Prong->GetEntriesFast();\r
- if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);\r
- \r
- \r
- Int_t pdgDstoKKpi[3]={321,321,211};\r
- Int_t nSelectedloose=0;\r
- Int_t nSelectedtight=0;\r
-\r
- for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {\r
- \r
- AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);\r
- fHistNEvents->Fill(8);\r
- \r
- if(fUseSelectionBit && !(d->HasSelectionBit(AliRDHFCuts::kDsCuts))){\r
- fHistNEvents->Fill(11);\r
- continue;\r
- }\r
- \r
- Bool_t unsetvtx=kFALSE;\r
- if(!d->GetOwnPrimaryVtx()){\r
- d->SetOwnPrimaryVtx(vtx1);\r
- unsetvtx=kTRUE;\r
- }\r
- \r
- Bool_t recVtx=kFALSE;\r
- AliAODVertex *origownvtx=0x0;\r
- Int_t retCodeProdCuts=fProdCuts->IsSelected(d,AliRDHFCuts::kCandidate,aod);\r
- \r
- if(retCodeProdCuts) {\r
- if(fProdCuts->GetIsPrimaryWithoutDaughters()){\r
- if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx()); \r
- if(fProdCuts->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;\r
- else fProdCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);\r
- }\r
- } \r
- \r
- Double_t ptCand = d->Pt();\r
- Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);\r
- Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kCandidate,aod);\r
- Double_t rapid=d->YDs(); \r
- fYVsPt->Fill(ptCand,rapid);\r
-\r
- Bool_t isFidAcc=fAnalysisCuts->IsInFiducialAcceptance(ptCand,rapid);\r
- \r
- if(retCodeProdCuts>0){\r
- if(isFidAcc){\r
- nSelectedloose++;\r
- fHistNEvents->Fill(9);\r
- if(retCodeAnalysisCuts>0)nSelectedtight++;\r
- }\r
- }\r
- \r
- if(retCodeAnalysisCuts<=0) continue;\r
- if(!isFidAcc) continue;\r
- fHistNEvents->Fill(10);\r
- \r
- Int_t index=GetHistoIndex(iPtBin);\r
- fPtCandHist[index]->Fill(ptCand);\r
-\r
- Int_t isKKpi=retCodeAnalysisCuts&1;\r
- Int_t ispiKK=retCodeAnalysisCuts&2;\r
- Int_t isPhiKKpi=retCodeAnalysisCuts&4;\r
- Int_t isPhipiKK=retCodeAnalysisCuts&8;\r
- Int_t isK0starKKpi=retCodeAnalysisCuts&16; \r
- Int_t isK0starpiKK=retCodeAnalysisCuts&32;\r
-\r
- fChanHist[0]->Fill(retCodeAnalysisCuts);\r
- \r
- Int_t indexMCKKpi=-1;\r
- Int_t indexMCpiKK=-1;\r
- Int_t labDs=-1;\r
- Int_t pdgCode0=-999;\r
- \r
- if(fReadMC){\r
- labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);\r
- if(labDs>=0){\r
- Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();\r
- AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));\r
- pdgCode0=TMath::Abs(p->GetPdgCode());\r
-\r
- if(isKKpi){\r
- if(pdgCode0==321) { \r
- indexMCKKpi=GetSignalHistoIndex(iPtBin);\r
- fYVsPtSig->Fill(ptCand,rapid);\r
- fChanHist[1]->Fill(retCodeAnalysisCuts);\r
- }else{\r
- indexMCKKpi=GetReflSignalHistoIndex(iPtBin);\r
- fChanHist[3]->Fill(retCodeAnalysisCuts);\r
- }\r
- }\r
- if(ispiKK){\r
- if(pdgCode0==211) { \r
- indexMCpiKK=GetSignalHistoIndex(iPtBin);\r
- fYVsPtSig->Fill(ptCand,rapid);\r
- fChanHist[1]->Fill(retCodeAnalysisCuts);\r
- }else{\r
- indexMCpiKK=GetReflSignalHistoIndex(iPtBin);\r
- fChanHist[3]->Fill(retCodeAnalysisCuts);\r
- }\r
- }\r
- }else{\r
- indexMCpiKK=GetBackgroundHistoIndex(iPtBin);\r
- indexMCKKpi=GetBackgroundHistoIndex(iPtBin);\r
- fChanHist[2]->Fill(retCodeAnalysisCuts);\r
- }\r
- }\r
-\r
- if(isKKpi){\r
- Double_t invMass=d->InvMassDsKKpi();\r
- fMassHist[index]->Fill(invMass);\r
- fPtVsMass->Fill(invMass,ptCand);\r
- if(isPhiKKpi) fMassHistPhi[index]->Fill(invMass); \r
- if(isK0starKKpi) fMassHistK0st[index]->Fill(invMass);\r
- if(fReadMC && indexMCKKpi!=-1){\r
- fMassHist[indexMCKKpi]->Fill(invMass);\r
- if(isPhiKKpi) fMassHistPhi[indexMCKKpi]->Fill(invMass);\r
- if(isK0starKKpi) fMassHistK0st[indexMCKKpi]->Fill(invMass); \r
- }\r
- }\r
- if(ispiKK){\r
- Double_t invMass=d->InvMassDspiKK();\r
- fMassHist[index]->Fill(invMass);\r
- fPtVsMass->Fill(invMass,ptCand);\r
- if(isPhipiKK) fMassHistPhi[index]->Fill(invMass);\r
- if(isK0starpiKK) fMassHistK0st[index]->Fill(invMass);\r
- if(fReadMC && indexMCpiKK!=-1){\r
- fMassHist[indexMCpiKK]->Fill(invMass);\r
- if(isPhipiKK) fMassHistPhi[indexMCpiKK]->Fill(invMass);\r
- if(isK0starpiKK) fMassHistK0st[indexMCpiKK]->Fill(invMass); \r
- }\r
- }\r
-\r
- if(fDoCutVarHistos){\r
- Double_t dlen=d->DecayLength();\r
- Double_t cosp=d->CosPointingAngle();\r
- Double_t pt0=d->PtProng(0);\r
- Double_t pt1=d->PtProng(1);\r
- Double_t pt2=d->PtProng(2);\r
- Double_t sigvert=d->GetSigmaVert();\r
- Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);\r
- Double_t dca=d->GetDCA(); \r
- Double_t ptmax=0;\r
- for(Int_t i=0;i<3;i++){\r
- if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);\r
- }\r
- fCosPHist[index]->Fill(cosp);\r
- fDLenHist[index]->Fill(dlen);\r
- fSigVertHist[index]->Fill(sigvert);\r
- fSumd02Hist[index]->Fill(sumD02);\r
- fPtMaxHist[index]->Fill(ptmax);\r
- fDCAHist[index]->Fill(dca);\r
- fPtProng0Hist[index]->Fill(pt0);\r
- fPtProng1Hist[index]->Fill(pt1);\r
- fPtProng2Hist[index]->Fill(pt2);\r
- if(isKKpi){\r
- Double_t massKK=d->InvMass2Prongs(0,1,321,321);\r
- Double_t massKp=d->InvMass2Prongs(1,2,321,211);\r
- fDalitz[index]->Fill(massKK,massKp);\r
- if(isPhiKKpi) fDalitzPhi[index]->Fill(massKK,massKp);\r
- if(isK0starKKpi) fDalitzK0st[index]->Fill(massKK,massKp);\r
- if(fReadMC && indexMCKKpi!=-1){\r
- fDalitz[indexMCKKpi]->Fill(massKK,massKp);\r
- if(isPhiKKpi) fDalitzPhi[indexMCKKpi]->Fill(massKK,massKp);\r
- if(isK0starKKpi) fDalitzK0st[indexMCKKpi]->Fill(massKK,massKp);\r
- fCosPHist[indexMCKKpi]->Fill(cosp);\r
- fDLenHist[indexMCKKpi]->Fill(dlen);\r
- fSigVertHist[indexMCKKpi]->Fill(sigvert);\r
- fSumd02Hist[indexMCKKpi]->Fill(sumD02);\r
- fPtMaxHist[indexMCKKpi]->Fill(ptmax);\r
- fPtCandHist[indexMCKKpi]->Fill(ptCand);\r
- fDCAHist[indexMCKKpi]->Fill(dca);\r
- fPtProng0Hist[indexMCKKpi]->Fill(pt0);\r
- fPtProng1Hist[indexMCKKpi]->Fill(pt1);\r
- fPtProng2Hist[indexMCKKpi]->Fill(pt2);\r
- } \r
- }\r
- if(ispiKK){\r
- Double_t massKK=d->InvMass2Prongs(1,2,321,321);\r
- Double_t massKp=d->InvMass2Prongs(0,1,211,321);\r
- fDalitz[index]->Fill(massKK,massKp);\r
- if(isPhipiKK) fDalitzPhi[index]->Fill(massKK,massKp);\r
- if(isK0starpiKK) fDalitzK0st[index]->Fill(massKK,massKp);\r
-\r
- \r
- if(fReadMC && indexMCpiKK!=-1){\r
- fDalitz[indexMCpiKK]->Fill(massKK,massKp);\r
- if(isPhipiKK) fDalitzPhi[indexMCpiKK]->Fill(massKK,massKp);\r
- if(isK0starpiKK) fDalitzK0st[indexMCpiKK]->Fill(massKK,massKp);\r
- fCosPHist[indexMCpiKK]->Fill(cosp);\r
- fDLenHist[indexMCpiKK]->Fill(dlen);\r
- fSigVertHist[indexMCpiKK]->Fill(sigvert);\r
- fSumd02Hist[indexMCpiKK]->Fill(sumD02);\r
- fPtMaxHist[indexMCpiKK]->Fill(ptmax);\r
- fPtCandHist[indexMCpiKK]->Fill(ptCand);\r
- fDCAHist[indexMCpiKK]->Fill(dca);\r
- fPtProng0Hist[indexMCpiKK]->Fill(pt0);\r
- fPtProng1Hist[indexMCpiKK]->Fill(pt1);\r
- fPtProng2Hist[indexMCpiKK]->Fill(pt2);\r
- }\r
- }\r
- \r
- \r
- }\r
- \r
- Float_t tmp[36];\r
- if(fFillNtuple>0){\r
- \r
- if ((fFillNtuple==1 && (isPhiKKpi || isPhipiKK)) || (fFillNtuple==2 && (isK0starKKpi || isK0starpiKK)) || (fFillNtuple==3 && (isKKpi || ispiKK))){\r
- \r
- AliAODTrack *track0=(AliAODTrack*)d->GetDaughter(0);\r
- AliAODTrack *track1=(AliAODTrack*)d->GetDaughter(1);\r
- AliAODTrack *track2=(AliAODTrack*)d->GetDaughter(2);\r
- \r
- UInt_t BitMapPIDTrack0=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track0);\r
- UInt_t BitMapPIDTrack1=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track1);\r
- UInt_t BitMapPIDTrack2=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track2);\r
- \r
- tmp[0]=Float_t(labDs);\r
- tmp[1]=Float_t(retCodeAnalysisCuts);\r
- tmp[2]=Float_t(pdgCode0); \r
- tmp[3]=d->PtProng(0);\r
- tmp[4]=d->PtProng(1);\r
- tmp[5]=d->PtProng(2);\r
- tmp[6]=d->Pt();\r
- tmp[7]=d->PProng(0);\r
- tmp[8]=d->PProng(1);\r
- tmp[9]=d->PProng(2);\r
- tmp[10]=Int_t(BitMapPIDTrack0);\r
- tmp[11]=Int_t(BitMapPIDTrack1);\r
- tmp[12]=Int_t(BitMapPIDTrack2);\r
- tmp[13]=d->CosPointingAngle();\r
- tmp[14]=d->CosPointingAngleXY();\r
- tmp[15]=d->DecayLength();\r
- tmp[16]=d->DecayLengthXY();\r
- tmp[17]=d->NormalizedDecayLength();\r
- tmp[18]=d->NormalizedDecayLengthXY();\r
- tmp[19]=d->InvMassDsKKpi();\r
- tmp[20]=d->InvMassDspiKK();\r
- tmp[21]=d->GetSigmaVert();\r
- tmp[22]=d->Getd0Prong(0);\r
- tmp[23]=d->Getd0Prong(1);\r
- tmp[24]=d->Getd0Prong(2);\r
- tmp[25]=d->GetDCA();\r
- tmp[26]=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);\r
- tmp[27]=d->InvMass2Prongs(0,1,321,321);\r
- tmp[28]=d->InvMass2Prongs(1,2,321,321);\r
- tmp[29]=d->InvMass2Prongs(1,2,321,211);\r
- tmp[30]=d->InvMass2Prongs(0,1,211,321);\r
- tmp[31]=d->CosPiDsLabFrameKKpi(); \r
- tmp[32]=d->CosPiDsLabFramepiKK(); \r
- tmp[33]=d->CosPiKPhiRFrameKKpi(); \r
- tmp[34]=d->CosPiKPhiRFramepiKK(); \r
- tmp[35]=(Float_t)(centrality);\r
- \r
- \r
- \r
- fNtupleDs->Fill(tmp);\r
- PostData(4,fNtupleDs);\r
- } \r
- }\r
- \r
- if(unsetvtx) d->UnsetOwnPrimaryVtx();\r
- if(recVtx)fProdCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);\r
- }\r
- \r
- fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);\r
- fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);\r
-\r
- PostData(1,fOutput); \r
- PostData(3,fCounter); \r
-\r
- return;\r
-}\r
-\r
-//_________________________________________________________________\r
-\r
-void AliAnalysisTaskSEDs::Terminate(Option_t */*option*/)\r
-{\r
- // Terminate analysis\r
- //\r
- if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");\r
- fOutput = dynamic_cast<TList*> (GetOutputData(1));\r
- if (!fOutput) { \r
- printf("ERROR: fOutput not available\n");\r
- return;\r
- }\r
- fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));\r
- if(fHistNEvents){\r
- printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));\r
- }else{\r
- printf("ERROR: fHistNEvents not available\n");\r
- return;\r
- }\r
- return;\r
-}\r
-\r
+/**************************************************************************
+ * Copyright(c) 2007-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 *
+ * appear 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$ */
+
+///////////////////////////////////////////////////////////////////
+// //
+// Analysis task to produce Ds candidates mass spectra //
+// Origin: F.Prino, Torino, prino@to.infn.it //
+// //
+///////////////////////////////////////////////////////////////////
+
+#include <TClonesArray.h>
+#include <TNtuple.h>
+#include <TList.h>
+#include <TString.h>
+#include <TH1F.h>
+#include <TMath.h>
+#include <TDatabasePDG.h>
+
+#include "AliAnalysisManager.h"
+#include "AliAODHandler.h"
+#include "AliAODEvent.h"
+#include "AliAODVertex.h"
+#include "AliAODTrack.h"
+#include "AliAODMCHeader.h"
+#include "AliAODMCParticle.h"
+#include "AliAODRecoDecayHF3Prong.h"
+#include "AliAnalysisVertexingHF.h"
+#include "AliRDHFCutsDstoKKpi.h"
+#include "AliAnalysisTaskSE.h"
+#include "AliNormalizationCounter.h"
+#include "AliAnalysisTaskSEDs.h"
+
+ClassImp(AliAnalysisTaskSEDs)
+
+
+//________________________________________________________________________
+AliAnalysisTaskSEDs::AliAnalysisTaskSEDs():
+ AliAnalysisTaskSE(),
+ fOutput(0),
+ fHistNEvents(0),
+ fPtVsMass(0),
+ fYVsPt(0),
+ fYVsPtSig(0),
+ fNtupleDs(0),
+ fFillNtuple(0),
+ fReadMC(kFALSE),
+ fDoCutVarHistos(kTRUE),
+ fUseSelectionBit(kFALSE),
+ fNPtBins(0),
+ fListCuts(0),
+ fMassRange(0.8),
+ fMassBinSize(0.002),
+ fCounter(0),
+ fAnalysisCuts(0)
+{
+ // Default constructor
+
+ for(Int_t i=0;i<4;i++) {
+ if(fChanHist[i]) fChanHist[i]=0;
+ }
+
+ for(Int_t i=0;i<4*kMaxPtBins;i++){
+
+ if(fPtCandHist[i]) fPtCandHist[i]=0;
+ if(fMassHist[i]) fMassHist[i]=0;
+ if(fCosPHist[i]) fCosPHist[i]=0;
+ if(fDLenHist[i]) fDLenHist[i]=0;
+ if(fSumd02Hist[i]) fSumd02Hist[i]=0;
+ if(fSigVertHist[i]) fSigVertHist[i]=0;
+ if(fPtMaxHist[i]) fPtMaxHist[i]=0;
+ if(fDCAHist[i]) fDCAHist[i]=0;
+ if(fPtProng0Hist[i]) fPtProng0Hist[i]=0;
+ if(fPtProng1Hist[i]) fPtProng1Hist[i]=0;
+ if(fPtProng2Hist[i]) fPtProng2Hist[i]=0;
+ if(fDalitz[i]) fDalitz[i]=0;
+ if(fDalitzPhi[i]) fDalitzPhi[i]=0;
+ if(fDalitzK0st[i]) fDalitzK0st[i]=0;
+
+ }
+ for(Int_t i=0;i<3*kMaxPtBins;i++){
+ if(fMassHistPhi[i]) fMassHistPhi[i]=0;
+ if(fMassHistK0st[i]) fMassHistK0st[i]=0;
+
+ }
+
+ for(Int_t i=0;i<kMaxPtBins+1;i++){
+ fPtLimits[i]=0;
+ }
+
+}
+
+//________________________________________________________________________
+AliAnalysisTaskSEDs::AliAnalysisTaskSEDs(const char *name,AliRDHFCutsDstoKKpi* analysiscuts,Int_t fillNtuple):
+ AliAnalysisTaskSE(name),
+ fOutput(0),
+ fHistNEvents(0),
+ fPtVsMass(0),
+ fYVsPt(0),
+ fYVsPtSig(0),
+ fNtupleDs(0),
+ fFillNtuple(fillNtuple),
+ fReadMC(kFALSE),
+ fDoCutVarHistos(kTRUE),
+ fUseSelectionBit(kFALSE),
+ fNPtBins(0),
+ fListCuts(0),
+ fMassRange(0.8),
+ fMassBinSize(0.002),
+ fCounter(0),
+ fAnalysisCuts(analysiscuts)
+{
+ // Default constructor
+ // Output slot #1 writes into a TList container
+
+ for(Int_t i=0;i<4;i++) {
+ if(fChanHist[i]) fChanHist[i]=0;
+ }
+
+ for(Int_t i=0;i<4*kMaxPtBins;i++){
+
+ if(fPtCandHist[i]) fPtCandHist[i]=0;
+ if(fMassHist[i]) fMassHist[i]=0;
+ if(fCosPHist[i]) fCosPHist[i]=0;
+ if(fDLenHist[i]) fDLenHist[i]=0;
+ if(fSumd02Hist[i]) fSumd02Hist[i]=0;
+ if(fSigVertHist[i]) fSigVertHist[i]=0;
+ if(fPtMaxHist[i]) fPtMaxHist[i]=0;
+ if(fDCAHist[i]) fDCAHist[i]=0;
+ if(fPtProng0Hist[i]) fPtProng0Hist[i]=0;
+ if(fPtProng1Hist[i]) fPtProng1Hist[i]=0;
+ if(fPtProng2Hist[i]) fPtProng2Hist[i]=0;
+ if(fDalitz[i]) fDalitz[i]=0;
+ if(fDalitzPhi[i]) fDalitzPhi[i]=0;
+ if(fDalitzK0st[i]) fDalitzK0st[i]=0;
+
+ }
+ for(Int_t i=0;i<3*kMaxPtBins;i++){
+ if(fMassHistPhi[i]) fMassHistPhi[i]=0;
+ if(fMassHistK0st[i]) fMassHistK0st[i]=0;
+
+ }
+
+ for(Int_t i=0;i<kMaxPtBins+1;i++){
+ fPtLimits[i]=0;
+ }
+
+ Int_t nptbins=fAnalysisCuts->GetNPtBins();
+ Float_t *ptlim=fAnalysisCuts->GetPtBinLimits();
+ SetPtBins(nptbins,ptlim);
+
+ DefineOutput(1,TList::Class()); //My private output
+
+ DefineOutput(2,TList::Class());
+
+ DefineOutput(3,AliNormalizationCounter::Class());
+
+ if(fFillNtuple>0){
+ // Output slot #4 writes into a TNtuple container
+ DefineOutput(4,TNtuple::Class()); //My private output
+ }
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSEDs::SetPtBins(Int_t n, Float_t* lim){
+ // define pt bins for analysis
+ if(n>kMaxPtBins){
+ printf("Max. number of Pt bins = %d\n",kMaxPtBins);
+ fNPtBins=kMaxPtBins;
+ fPtLimits[0]=0.;
+ fPtLimits[1]=1.;
+ fPtLimits[2]=3.;
+ fPtLimits[3]=5.;
+ fPtLimits[4]=10.;
+ for(Int_t i=5; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
+ }else{
+ fNPtBins=n;
+ for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
+ for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
+ }
+ if(fDebug > 1){
+ printf("Number of Pt bins = %d\n",fNPtBins);
+ for(Int_t i=0; i<fNPtBins; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fPtLimits[i],fPtLimits[i+1]);
+ }
+}
+//________________________________________________________________________
+AliAnalysisTaskSEDs::~AliAnalysisTaskSEDs()
+{
+ // Destructor
+ delete fOutput;
+ delete fHistNEvents;
+ delete fListCuts;
+
+ for(Int_t i=0;i<4*fNPtBins;i++){
+
+ if(fMassHist[i]){ delete fMassHist[i]; fMassHist[i]=0;}
+ if(fMassHistPhi[i]){ delete fMassHistPhi[i]; fMassHistPhi[i]=0;}
+ if(fMassHistK0st[i]){ delete fMassHistK0st[i]; fMassHistK0st[i]=0;}
+ if(fCosPHist[i]){ delete fCosPHist[i]; fCosPHist[i]=0;}
+ if(fDLenHist[i]){ delete fDLenHist[i]; fDLenHist[i]=0;}
+ if(fSumd02Hist[i]){ delete fSumd02Hist[i]; fSumd02Hist[i]=0;}
+ if(fSigVertHist[i]){ delete fSigVertHist[i]; fSigVertHist[i]=0;}
+ if(fPtMaxHist[i]){ delete fPtMaxHist[i]; fPtMaxHist[i]=0;}
+ if(fDCAHist[i]){ delete fDCAHist[i]; fDCAHist[i]=0;}
+ if(fPtProng0Hist[i]){ delete fPtProng0Hist[i]; fPtProng0Hist[i]=0;}
+ if(fPtProng1Hist[i]){ delete fPtProng1Hist[i]; fPtProng1Hist[i]=0;}
+ if(fPtProng2Hist[i]){ delete fPtProng2Hist[i]; fPtProng2Hist[i]=0;}
+ if(fDalitz[i]){ delete fDalitz[i]; fDalitz[i]=0;}
+ if(fDalitzPhi[i]){ delete fDalitzPhi[i]; fDalitzPhi[i]=0;}
+ if(fDalitzK0st[i]){ delete fDalitzK0st[i]; fDalitzK0st[i]=0;}
+
+ }
+
+ delete fPtVsMass;
+ delete fYVsPt;
+ delete fYVsPtSig;
+ delete fNtupleDs;
+ delete fCounter;
+ delete fAnalysisCuts;
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSEDs::Init()
+{
+ // Initialization
+
+ if(fDebug > 1) printf("AnalysisTaskSEDs::Init() \n");
+
+ fListCuts=new TList();
+ fListCuts->SetOwner();
+ fListCuts->SetName("CutObjects");
+
+ AliRDHFCutsDstoKKpi *analysis = new AliRDHFCutsDstoKKpi(*fAnalysisCuts);
+ analysis->SetName("AnalysisCuts");
+
+ fListCuts->Add(analysis);
+ PostData(2,fListCuts);
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSEDs::UserCreateOutputObjects()
+{
+ // Create the output container
+ //
+ if(fDebug > 1) printf("AnalysisTaskSEDs::UserCreateOutputObjects() \n");
+
+ // Several histograms are more conveniently managed in a TList
+ fOutput = new TList();
+ fOutput->SetOwner();
+ fOutput->SetName("OutputHistos");
+
+ Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
+
+ Int_t nInvMassBins=(Int_t)(fMassRange/fMassBinSize+0.5);
+ if(nInvMassBins%2==1) nInvMassBins++;
+ // Double_t minMass=massDs-1.0*nInvMassBins*fMassBinSize;
+ Double_t minMass=massDs-0.5*nInvMassBins*fMassBinSize;
+ // Double_t maxMass=massDs+1.0*nInvMassBins*fMassBinSize;
+ Double_t maxMass=massDs+0.5*nInvMassBins*fMassBinSize;
+
+ TString hisname;
+ TString htype;
+ Int_t index;
+ for(Int_t iType=0; iType<4; iType++){
+ for(Int_t i=0;i<fNPtBins;i++){
+ if(iType==0){
+ htype="All";
+ index=GetHistoIndex(i);
+ }else if(iType==1){
+ htype="Sig";
+ index=GetSignalHistoIndex(i);
+ }else if(iType==2){
+ htype="Bkg";
+ index=GetBackgroundHistoIndex(i);
+ }else{
+ htype="ReflSig";
+ index=GetReflSignalHistoIndex(i);
+ }
+ hisname.Form("hMass%sPt%d",htype.Data(),i);
+ fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
+ fMassHist[index]->Sumw2();
+ hisname.Form("hMass%sPt%dphi",htype.Data(),i);
+ fMassHistPhi[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
+ fMassHistPhi[index]->Sumw2();
+ hisname.Form("hMass%sPt%dk0st",htype.Data(),i);
+ fMassHistK0st[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
+ fMassHistK0st[index]->Sumw2();
+ hisname.Form("hCosP%sPt%d",htype.Data(),i);
+ fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
+ fCosPHist[index]->Sumw2();
+ hisname.Form("hDLen%sPt%d",htype.Data(),i);
+ fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
+ fDLenHist[index]->Sumw2();
+ hisname.Form("hSumd02%sPt%d",htype.Data(),i);
+ fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
+ fSumd02Hist[index]->Sumw2();
+ hisname.Form("hSigVert%sPt%d",htype.Data(),i);
+ fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
+ fSigVertHist[index]->Sumw2();
+ hisname.Form("hPtMax%sPt%d",htype.Data(),i);
+ fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
+ fPtMaxHist[index]->Sumw2();
+ hisname.Form("hPtCand%sPt%d",htype.Data(),i);
+ fPtCandHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
+ fPtCandHist[index]->Sumw2();
+ hisname.Form("hDCA%sPt%d",htype.Data(),i);
+ fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
+ fDCAHist[index]->Sumw2();
+ hisname.Form("hPtProng0%sPt%d",htype.Data(),i);
+ fPtProng0Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
+ fPtProng0Hist[index]->Sumw2();
+ hisname.Form("hPtProng1%sPt%d",htype.Data(),i);
+ fPtProng1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
+ fPtProng1Hist[index]->Sumw2();
+ hisname.Form("hPtProng2%sPt%d",htype.Data(),i);
+ fPtProng2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
+ fPtProng2Hist[index]->Sumw2();
+ hisname.Form("hDalitz%sPt%d",htype.Data(),i);
+ fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
+ fDalitz[index]->Sumw2();
+ hisname.Form("hDalitz%sPt%dphi",htype.Data(),i);
+ fDalitzPhi[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
+ fDalitzPhi[index]->Sumw2();
+ hisname.Form("hDalitz%sPt%dk0st",htype.Data(),i);
+ fDalitzK0st[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
+ fDalitzK0st[index]->Sumw2();
+ }
+ }
+
+ for(Int_t i=0; i<4*fNPtBins; i++){
+ fOutput->Add(fMassHist[i]);
+ fOutput->Add(fMassHistPhi[i]);
+ fOutput->Add(fMassHistK0st[i]);
+ fOutput->Add(fCosPHist[i]);
+ fOutput->Add(fDLenHist[i]);
+ fOutput->Add(fSumd02Hist[i]);
+ fOutput->Add(fSigVertHist[i]);
+ fOutput->Add(fPtMaxHist[i]);
+ fOutput->Add(fPtCandHist[i]);
+ fOutput->Add(fDCAHist[i]);
+ fOutput->Add(fPtProng0Hist[i]);
+ fOutput->Add(fPtProng1Hist[i]);
+ fOutput->Add(fPtProng2Hist[i]);
+ fOutput->Add(fDalitz[i]);
+ fOutput->Add(fDalitzPhi[i]);
+ fOutput->Add(fDalitzK0st[i]);
+ }
+
+ fChanHist[0] = new TH1F("hChanAll", "KKpi and piKK candidates",64,-0.5,63.5);
+ fChanHist[1] = new TH1F("hChanSig", "KKpi and piKK candidates",64,-0.5,63.5);
+ fChanHist[2] = new TH1F("hChanBkg", "KKpi and piKK candidates",64,-0.5,63.5);
+ fChanHist[3] = new TH1F("hChanReflSig", "KKpi and piKK candidates",64,-0.5,63.5);
+ for(Int_t i=0;i<4;i++){
+ fChanHist[i]->Sumw2();
+ fChanHist[i]->SetMinimum(0);
+ fOutput->Add(fChanHist[i]);
+ }
+
+ fHistNEvents = new TH1F("hNEvents", "number of events ",11,-0.5,10.5);
+ fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
+ fHistNEvents->GetXaxis()->SetBinLabel(2,"n. passing IsEvSelected");
+ fHistNEvents->GetXaxis()->SetBinLabel(3,"n. rejected due to trigger");
+ fHistNEvents->GetXaxis()->SetBinLabel(4,"n. rejected due to not reco vertex");
+ fHistNEvents->GetXaxis()->SetBinLabel(5,"n. rejected for contr vertex");
+ fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected for vertex out of accept");
+ fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected for pileup events");
+ fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events");
+ fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of 3 prong candidates");
+ fHistNEvents->GetXaxis()->SetBinLabel(10,"no. of Ds after filtering cuts");
+ fHistNEvents->GetXaxis()->SetBinLabel(11,"no. of Ds after selection cuts");
+
+ fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
+
+ fHistNEvents->Sumw2();
+ fHistNEvents->SetMinimum(0);
+ fOutput->Add(fHistNEvents);
+
+ fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nInvMassBins,minMass,maxMass,40,0.,20.);
+ fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
+ fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
+
+ fOutput->Add(fPtVsMass);
+ fOutput->Add(fYVsPt);
+ fOutput->Add(fYVsPtSig);
+
+ //Counter for Normalization
+ fCounter = new AliNormalizationCounter("NormalizationCounter");
+ fCounter->Init();
+
+ PostData(1,fOutput);
+ PostData(3,fCounter);
+
+ if(fFillNtuple>0){
+ OpenFile(4); // 4 is the slot number of the ntuple
+
+ fNtupleDs = new TNtuple("fNtupleDs","Ds","labDs:retcode:pdgcode0:Pt0:Pt1:Pt2:PtRec:P0:P1:P2:PidTrackBit0:PidTrackBit1:PidTrackBit2:PointingAngle:PointingAngleXY:DecLeng:DecLengXY:NorDecLeng:NorDecLengXY:InvMassKKpi:InvMasspiKK:sigvert:d00:d01:d02:dca:d0square:InvMassPhiKKpi:InvMassPhipiKK:InvMassK0starKKpi:InvMassK0starpiKK:cosinePiDsFrameKKpi:cosinePiDsFramepiKK:cosineKPhiFrameKKpi:cosineKPhiFramepiKK:centrality");
+
+ }
+
+
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSEDs::UserExec(Option_t */*option*/)
+{
+ // Ds selection for current event, fill mass histos and selecetion variable histo
+ // separate signal and backgound if fReadMC is activated
+
+ AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
+
+ TClonesArray *array3Prong = 0;
+ if(!aod && 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.
+ aod = 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();
+ array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
+ }
+ } else if(aod) {
+ array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
+ }
+
+ if(!aod || !array3Prong) {
+ printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
+ return;
+ }
+
+
+ // fix for temporary bug in ESDfilter
+ // the AODs with null vertex pointer didn't pass the PhysSel
+ if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
+
+
+ fHistNEvents->Fill(0); // count event
+ // Post the data already here
+ PostData(1,fOutput);
+
+ fCounter->StoreEvent(aod,fAnalysisCuts,fReadMC);
+
+
+ Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
+ if(fAnalysisCuts->IsEventRejectedDueToTrigger())fHistNEvents->Fill(2);
+ if(fAnalysisCuts->IsEventRejectedDueToNotRecoVertex())fHistNEvents->Fill(3);
+ if(fAnalysisCuts->IsEventRejectedDueToVertexContributors())fHistNEvents->Fill(4);
+ if(fAnalysisCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion())fHistNEvents->Fill(5);
+ if(fAnalysisCuts->IsEventRejectedDueToPileupSPD())fHistNEvents->Fill(6);
+ if(fAnalysisCuts->IsEventRejectedDueToCentrality())fHistNEvents->Fill(7);
+
+ Float_t centrality=fAnalysisCuts->GetCentrality(aod);
+ //Int_t runNumber=aod->GetRunNumber();
+
+ if(!isEvSel)return;
+
+
+
+
+ fHistNEvents->Fill(1);
+
+ TClonesArray *arrayMC=0;
+ AliAODMCHeader *mcHeader=0;
+
+ // AOD primary vertex
+ AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
+ // vtx1->Print();
+
+ // load MC particles
+ if(fReadMC){
+
+ arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
+ if(!arrayMC) {
+ printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
+ return;
+ }
+
+ // load MC header
+ mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
+ if(!mcHeader) {
+ printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
+ return;
+ }
+ }
+
+ Int_t n3Prong = array3Prong->GetEntriesFast();
+ if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
+
+
+ Int_t pdgDstoKKpi[3]={321,321,211};
+ Int_t nSelected=0;
+ Int_t nFiltered=0;
+ for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
+
+ AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
+ fHistNEvents->Fill(8);
+
+ if(fUseSelectionBit && !(d->HasSelectionBit(AliRDHFCuts::kDsCuts))){
+ continue;
+ }
+ nFiltered++;
+ fHistNEvents->Fill(9);
+
+ Bool_t unsetvtx=kFALSE;
+ if(!d->GetOwnPrimaryVtx()){
+ d->SetOwnPrimaryVtx(vtx1);
+ unsetvtx=kTRUE;
+ }
+
+ Bool_t recVtx=kFALSE;
+ AliAODVertex *origownvtx=0x0;
+
+
+ Double_t ptCand = d->Pt();
+ Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);
+ Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kCandidate,aod);
+ Double_t rapid=d->YDs();
+ fYVsPt->Fill(ptCand,rapid);
+
+ if(retCodeAnalysisCuts<=0) continue;
+ Bool_t isFidAcc=fAnalysisCuts->IsInFiducialAcceptance(ptCand,rapid);
+ if(!isFidAcc) continue;
+
+ if(fAnalysisCuts->GetIsPrimaryWithoutDaughters()){
+ if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
+ if(fAnalysisCuts->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
+ else fAnalysisCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
+ }
+
+
+ fHistNEvents->Fill(10);
+ nSelected++;
+
+ Int_t index=GetHistoIndex(iPtBin);
+ fPtCandHist[index]->Fill(ptCand);
+
+ Int_t isKKpi=retCodeAnalysisCuts&1;
+ Int_t ispiKK=retCodeAnalysisCuts&2;
+ Int_t isPhiKKpi=retCodeAnalysisCuts&4;
+ Int_t isPhipiKK=retCodeAnalysisCuts&8;
+ Int_t isK0starKKpi=retCodeAnalysisCuts&16;
+ Int_t isK0starpiKK=retCodeAnalysisCuts&32;
+
+ fChanHist[0]->Fill(retCodeAnalysisCuts);
+
+ Int_t indexMCKKpi=-1;
+ Int_t indexMCpiKK=-1;
+ Int_t labDs=-1;
+ Int_t pdgCode0=-999;
+
+ if(fReadMC){
+ labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
+ if(labDs>=0){
+ Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
+ AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));
+ pdgCode0=TMath::Abs(p->GetPdgCode());
+
+ if(isKKpi){
+ if(pdgCode0==321) {
+ indexMCKKpi=GetSignalHistoIndex(iPtBin);
+ fYVsPtSig->Fill(ptCand,rapid);
+ fChanHist[1]->Fill(retCodeAnalysisCuts);
+ }else{
+ indexMCKKpi=GetReflSignalHistoIndex(iPtBin);
+ fChanHist[3]->Fill(retCodeAnalysisCuts);
+ }
+ }
+ if(ispiKK){
+ if(pdgCode0==211) {
+ indexMCpiKK=GetSignalHistoIndex(iPtBin);
+ fYVsPtSig->Fill(ptCand,rapid);
+ fChanHist[1]->Fill(retCodeAnalysisCuts);
+ }else{
+ indexMCpiKK=GetReflSignalHistoIndex(iPtBin);
+ fChanHist[3]->Fill(retCodeAnalysisCuts);
+ }
+ }
+ }else{
+ indexMCpiKK=GetBackgroundHistoIndex(iPtBin);
+ indexMCKKpi=GetBackgroundHistoIndex(iPtBin);
+ fChanHist[2]->Fill(retCodeAnalysisCuts);
+ }
+ }
+
+ if(isKKpi){
+ Double_t invMass=d->InvMassDsKKpi();
+ fMassHist[index]->Fill(invMass);
+ fPtVsMass->Fill(invMass,ptCand);
+ if(isPhiKKpi) fMassHistPhi[index]->Fill(invMass);
+ if(isK0starKKpi) fMassHistK0st[index]->Fill(invMass);
+ if(fReadMC && indexMCKKpi!=-1){
+ fMassHist[indexMCKKpi]->Fill(invMass);
+ if(isPhiKKpi) fMassHistPhi[indexMCKKpi]->Fill(invMass);
+ if(isK0starKKpi) fMassHistK0st[indexMCKKpi]->Fill(invMass);
+ }
+ }
+ if(ispiKK){
+ Double_t invMass=d->InvMassDspiKK();
+ fMassHist[index]->Fill(invMass);
+ fPtVsMass->Fill(invMass,ptCand);
+ if(isPhipiKK) fMassHistPhi[index]->Fill(invMass);
+ if(isK0starpiKK) fMassHistK0st[index]->Fill(invMass);
+ if(fReadMC && indexMCpiKK!=-1){
+ fMassHist[indexMCpiKK]->Fill(invMass);
+ if(isPhipiKK) fMassHistPhi[indexMCpiKK]->Fill(invMass);
+ if(isK0starpiKK) fMassHistK0st[indexMCpiKK]->Fill(invMass);
+ }
+ }
+
+ if(fDoCutVarHistos){
+ Double_t dlen=d->DecayLength();
+ Double_t cosp=d->CosPointingAngle();
+ Double_t pt0=d->PtProng(0);
+ Double_t pt1=d->PtProng(1);
+ Double_t pt2=d->PtProng(2);
+ Double_t sigvert=d->GetSigmaVert();
+ Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
+ Double_t dca=d->GetDCA();
+ Double_t ptmax=0;
+ for(Int_t i=0;i<3;i++){
+ if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
+ }
+ fCosPHist[index]->Fill(cosp);
+ fDLenHist[index]->Fill(dlen);
+ fSigVertHist[index]->Fill(sigvert);
+ fSumd02Hist[index]->Fill(sumD02);
+ fPtMaxHist[index]->Fill(ptmax);
+ fDCAHist[index]->Fill(dca);
+ fPtProng0Hist[index]->Fill(pt0);
+ fPtProng1Hist[index]->Fill(pt1);
+ fPtProng2Hist[index]->Fill(pt2);
+ if(isKKpi){
+ Double_t massKK=d->InvMass2Prongs(0,1,321,321);
+ Double_t massKp=d->InvMass2Prongs(1,2,321,211);
+ fDalitz[index]->Fill(massKK,massKp);
+ if(isPhiKKpi) fDalitzPhi[index]->Fill(massKK,massKp);
+ if(isK0starKKpi) fDalitzK0st[index]->Fill(massKK,massKp);
+ if(fReadMC && indexMCKKpi!=-1){
+ fDalitz[indexMCKKpi]->Fill(massKK,massKp);
+ if(isPhiKKpi) fDalitzPhi[indexMCKKpi]->Fill(massKK,massKp);
+ if(isK0starKKpi) fDalitzK0st[indexMCKKpi]->Fill(massKK,massKp);
+ fCosPHist[indexMCKKpi]->Fill(cosp);
+ fDLenHist[indexMCKKpi]->Fill(dlen);
+ fSigVertHist[indexMCKKpi]->Fill(sigvert);
+ fSumd02Hist[indexMCKKpi]->Fill(sumD02);
+ fPtMaxHist[indexMCKKpi]->Fill(ptmax);
+ fPtCandHist[indexMCKKpi]->Fill(ptCand);
+ fDCAHist[indexMCKKpi]->Fill(dca);
+ fPtProng0Hist[indexMCKKpi]->Fill(pt0);
+ fPtProng1Hist[indexMCKKpi]->Fill(pt1);
+ fPtProng2Hist[indexMCKKpi]->Fill(pt2);
+ }
+ }
+ if(ispiKK){
+ Double_t massKK=d->InvMass2Prongs(1,2,321,321);
+ Double_t massKp=d->InvMass2Prongs(0,1,211,321);
+ fDalitz[index]->Fill(massKK,massKp);
+ if(isPhipiKK) fDalitzPhi[index]->Fill(massKK,massKp);
+ if(isK0starpiKK) fDalitzK0st[index]->Fill(massKK,massKp);
+
+
+ if(fReadMC && indexMCpiKK!=-1){
+ fDalitz[indexMCpiKK]->Fill(massKK,massKp);
+ if(isPhipiKK) fDalitzPhi[indexMCpiKK]->Fill(massKK,massKp);
+ if(isK0starpiKK) fDalitzK0st[indexMCpiKK]->Fill(massKK,massKp);
+ fCosPHist[indexMCpiKK]->Fill(cosp);
+ fDLenHist[indexMCpiKK]->Fill(dlen);
+ fSigVertHist[indexMCpiKK]->Fill(sigvert);
+ fSumd02Hist[indexMCpiKK]->Fill(sumD02);
+ fPtMaxHist[indexMCpiKK]->Fill(ptmax);
+ fPtCandHist[indexMCpiKK]->Fill(ptCand);
+ fDCAHist[indexMCpiKK]->Fill(dca);
+ fPtProng0Hist[indexMCpiKK]->Fill(pt0);
+ fPtProng1Hist[indexMCpiKK]->Fill(pt1);
+ fPtProng2Hist[indexMCpiKK]->Fill(pt2);
+ }
+ }
+
+
+ }
+
+ Float_t tmp[36];
+ if(fFillNtuple>0){
+
+ if ((fFillNtuple==1 && (isPhiKKpi || isPhipiKK)) || (fFillNtuple==2 && (isK0starKKpi || isK0starpiKK)) || (fFillNtuple==3 && (isKKpi || ispiKK))){
+
+ AliAODTrack *track0=(AliAODTrack*)d->GetDaughter(0);
+ AliAODTrack *track1=(AliAODTrack*)d->GetDaughter(1);
+ AliAODTrack *track2=(AliAODTrack*)d->GetDaughter(2);
+
+ UInt_t BitMapPIDTrack0=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track0);
+ UInt_t BitMapPIDTrack1=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track1);
+ UInt_t BitMapPIDTrack2=fAnalysisCuts->GetPIDTrackTPCTOFBitMap(track2);
+
+ tmp[0]=Float_t(labDs);
+ tmp[1]=Float_t(retCodeAnalysisCuts);
+ tmp[2]=Float_t(pdgCode0);
+ tmp[3]=d->PtProng(0);
+ tmp[4]=d->PtProng(1);
+ tmp[5]=d->PtProng(2);
+ tmp[6]=d->Pt();
+ tmp[7]=d->PProng(0);
+ tmp[8]=d->PProng(1);
+ tmp[9]=d->PProng(2);
+ tmp[10]=Int_t(BitMapPIDTrack0);
+ tmp[11]=Int_t(BitMapPIDTrack1);
+ tmp[12]=Int_t(BitMapPIDTrack2);
+ tmp[13]=d->CosPointingAngle();
+ tmp[14]=d->CosPointingAngleXY();
+ tmp[15]=d->DecayLength();
+ tmp[16]=d->DecayLengthXY();
+ tmp[17]=d->NormalizedDecayLength();
+ tmp[18]=d->NormalizedDecayLengthXY();
+ tmp[19]=d->InvMassDsKKpi();
+ tmp[20]=d->InvMassDspiKK();
+ tmp[21]=d->GetSigmaVert();
+ tmp[22]=d->Getd0Prong(0);
+ tmp[23]=d->Getd0Prong(1);
+ tmp[24]=d->Getd0Prong(2);
+ tmp[25]=d->GetDCA();
+ tmp[26]=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
+ tmp[27]=d->InvMass2Prongs(0,1,321,321);
+ tmp[28]=d->InvMass2Prongs(1,2,321,321);
+ tmp[29]=d->InvMass2Prongs(1,2,321,211);
+ tmp[30]=d->InvMass2Prongs(0,1,211,321);
+ tmp[31]=d->CosPiDsLabFrameKKpi();
+ tmp[32]=d->CosPiDsLabFramepiKK();
+ tmp[33]=d->CosPiKPhiRFrameKKpi();
+ tmp[34]=d->CosPiKPhiRFramepiKK();
+ tmp[35]=(Float_t)(centrality);
+
+
+
+ fNtupleDs->Fill(tmp);
+ PostData(4,fNtupleDs);
+ }
+ }
+
+ if(unsetvtx) d->UnsetOwnPrimaryVtx();
+ if(recVtx)fAnalysisCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
+ }
+
+ fCounter->StoreCandidates(aod,nFiltered,kTRUE);
+ fCounter->StoreCandidates(aod,nSelected,kFALSE);
+
+ PostData(1,fOutput);
+ PostData(3,fCounter);
+
+ return;
+}
+
+//_________________________________________________________________
+
+void AliAnalysisTaskSEDs::Terminate(Option_t */*option*/)
+{
+ // Terminate analysis
+ //
+ if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
+ fOutput = dynamic_cast<TList*> (GetOutputData(1));
+ if (!fOutput) {
+ printf("ERROR: fOutput not available\n");
+ return;
+ }
+ fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
+ if(fHistNEvents){
+ printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
+ }else{
+ printf("ERROR: fHistNEvents not available\n");
+ return;
+ }
+ return;
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff