--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-2008, 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. *
+ **************************************************************************/
+
+/////////////////////////////////////////////////////////////
+//
+// AliAnalysisTaskSE for the comparison of heavy flavor
+// decay candidates with the MC truth.
+//
+// Author: A.Dainese, andrea.dainese@lnl.infn.it
+/////////////////////////////////////////////////////////////
+
+#include <TClonesArray.h>
+#include <TNtuple.h>
+#include <TList.h>
+#include <TH1F.h>
+
+#include "AliAODEvent.h"
+#include "AliAODVertex.h"
+#include "AliAODTrack.h"
+#include "AliAODMCHeader.h"
+#include "AliAODMCParticle.h"
+#include "AliAODRecoDecayHF2Prong.h"
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisTaskSECompareHF.h"
+
+ClassImp(AliAnalysisTaskSECompareHF)
+
+
+//________________________________________________________________________
+AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF():
+AliAnalysisTaskSE(),
+fOutput(0),
+fNtupleD0Cmp(0),
+fHistMass(0)
+{
+ // Default constructor
+ SetD0toKpiCuts();
+
+ // Output slot #1 writes into a TList container
+ DefineOutput(1,TList::Class()); //My private output
+}
+
+//________________________________________________________________________
+AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name):
+AliAnalysisTaskSE(name),
+fOutput(0),
+fNtupleD0Cmp(0),
+fHistMass(0)
+{
+ // Default constructor
+ SetD0toKpiCuts();
+
+ // Output slot #1 writes into a TList container
+ DefineOutput(1,TList::Class()); //My private output
+}
+
+//________________________________________________________________________
+AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF()
+{
+ // Destructor
+ if (fOutput) {
+ delete fOutput;
+ fOutput = 0;
+ }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSECompareHF::Init()
+{
+ // Initialization
+
+ if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n");
+
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSECompareHF::UserCreateOutputObjects()
+{
+ // Create the output container
+ //
+ if(fDebug > 1) printf("AnalysisTaskSECompareHF::UserCreateOutputObjects() \n");
+
+ // Several histograms are more conveniently managed in a TList
+ fOutput = new TList();
+ fOutput->SetOwner();
+
+ fHistMass = new TH1F("fHistMass", "D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
+ fHistMass->Sumw2();
+ fHistMass->SetMinimum(0);
+ fOutput->Add(fHistMass);
+
+ fNtupleD0Cmp = new TNtuple("fNtupleD0Cmp","D0 comparison","pdg:VxRec:VxTrue:PtRec:PtTrue");
+ fOutput->Add(fNtupleD0Cmp);
+
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
+{
+ // Execute analysis for current event:
+ // heavy flavor candidates association to MC truth
+
+ AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
+
+ // load D0->Kpi candidates
+ TClonesArray *inputArrayD0toKpi =
+ (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
+ if(!inputArrayD0toKpi) {
+ printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
+ return;
+ }
+
+ // AOD primary vertex
+ AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
+ //vtx1->Print();
+
+ // load MC particles
+ TClonesArray *mcArray =
+ (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
+ if(!mcArray) {
+ printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n");
+ return;
+ }
+
+ // load MC header
+ AliAODMCHeader *mcHeader =
+ (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
+ if(!mcHeader) {
+ printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n");
+ return;
+ }
+
+
+
+ // loop over D0->Kpi candidates
+ Int_t nInD0toKpi = inputArrayD0toKpi->GetEntriesFast();
+ printf("Number of D0->Kpi: %d\n",nInD0toKpi);
+
+ Int_t lab0,lab1,labMother,labD0daugh0,labD0daugh1,pdgMother,pdgD0;
+
+ for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
+ AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayD0toKpi->UncheckedAt(iD0toKpi);
+ labD0daugh0=-1;
+ labD0daugh1=-1;
+ Bool_t unsetvtx=kFALSE;
+ if(!d->GetOwnPrimaryVtx()) {
+ d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
+ unsetvtx=kTRUE;
+ }
+ Int_t okD0=0,okD0bar=0;
+ if(d->SelectD0(fD0toKpiCuts,okD0,okD0bar)) {
+ // get daughter AOD tracks
+ AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
+ AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1);
+
+ lab0 = trk0->GetLabel();
+ lab1 = trk1->GetLabel();
+
+
+ AliAODMCParticle *part0 = (AliAODMCParticle*)mcArray->At(lab0);
+ if(!part0) {
+ printf("no MC particle\n");
+ continue;
+ }
+ while(part0->GetMother()>=0) {
+ labMother=part0->GetMother();
+ part0 = (AliAODMCParticle*)mcArray->At(labMother);
+ if(!part0) {
+ printf("no MC mother particle\n");
+ break;
+ }
+ pdgMother = TMath::Abs(part0->GetPdgCode());
+ if(pdgMother==421) {
+ labD0daugh0=labMother;
+ break;
+ }
+ }
+
+ AliAODMCParticle *part1 = (AliAODMCParticle*)mcArray->At(lab1);
+ if(!part1) {
+ printf("no MC particle\n");
+ continue;
+ }
+ while(part1->GetMother()>=0) {
+ labMother=part1->GetMother();
+ part1 = (AliAODMCParticle*)mcArray->At(labMother);
+ if(!part1) {
+ printf("no MC mother particle\n");
+ break;
+ }
+ pdgMother = TMath::Abs(part1->GetPdgCode());
+ if(pdgMother==421) {
+ labD0daugh1=labMother;
+ break;
+ }
+ }
+
+ // check if the candidate is signal
+ if(labD0daugh0>=0 && labD0daugh1>=0 && labD0daugh0==labD0daugh1) {
+
+ AliAODMCParticle *partD0 = (AliAODMCParticle*)mcArray->At(labD0daugh0);
+ pdgD0 = partD0->GetPdgCode();
+ Double_t invmass = (pdgD0==421 ? d->InvMassD0() : d->InvMassD0bar());
+ fHistMass->Fill(invmass);
+
+ // Post the data already here
+ PostData(1,fOutput);
+
+ fNtupleD0Cmp->Fill(pdgD0,d->Xv(),partD0->Xv(),d->Pt(),partD0->Pt());
+ }
+
+
+ }
+ if(unsetvtx) d->UnsetOwnPrimaryVtx();
+ } // end loop on D0->Kpi
+
+
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSECompareHF::Terminate(Option_t */*option*/)
+{
+ // Terminate analysis
+ //
+ if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n");
+
+ fOutput = dynamic_cast<TList*> (GetOutputData(1));
+ if (!fOutput) {
+ printf("ERROR: fOutput not available\n");
+ return;
+ }
+
+ fNtupleD0Cmp = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleD0Cmp"));
+ fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
+
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSECompareHF::SetD0toKpiCuts(Double_t cut0,Double_t cut1,
+ Double_t cut2,Double_t cut3,Double_t cut4,
+ Double_t cut5,Double_t cut6,
+ Double_t cut7,Double_t cut8)
+{
+ // Set the cuts for D0 selection
+ // cuts[0] = inv. mass half width [GeV]
+ // cuts[1] = dca [cm]
+ // cuts[2] = cosThetaStar
+ // cuts[3] = pTK [GeV/c]
+ // cuts[4] = pTPi [GeV/c]
+ // cuts[5] = d0K [cm] upper limit!
+ // cuts[6] = d0Pi [cm] upper limit!
+ // cuts[7] = d0d0 [cm^2]
+ // cuts[8] = cosThetaPoint
+
+ fD0toKpiCuts[0] = cut0;
+ fD0toKpiCuts[1] = cut1;
+ fD0toKpiCuts[2] = cut2;
+ fD0toKpiCuts[3] = cut3;
+ fD0toKpiCuts[4] = cut4;
+ fD0toKpiCuts[5] = cut5;
+ fD0toKpiCuts[6] = cut6;
+ fD0toKpiCuts[7] = cut7;
+ fD0toKpiCuts[8] = cut8;
+
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSECompareHF::SetD0toKpiCuts(const Double_t cuts[9])
+{
+ // Set the cuts for D0 selection
+ // cuts[0] = inv. mass half width [GeV]
+ // cuts[1] = dca [cm]
+ // cuts[2] = cosThetaStar
+ // cuts[3] = pTK [GeV/c]
+ // cuts[4] = pTPi [GeV/c]
+ // cuts[5] = d0K [cm] upper limit!
+ // cuts[6] = d0Pi [cm] upper limit!
+ // cuts[7] = d0d0 [cm^2]
+ // cuts[8] = cosThetaPoint
+
+ for(Int_t i=0; i<9; i++) fD0toKpiCuts[i] = cuts[i];
+
+ return;
+}
--- /dev/null
+#ifndef ALIANALYSISTASKCOMPAREHF_H
+#define ALIANALYSISTASKCOMPAREHF_H
+
+/* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+//*************************************************************************
+// Class AliAnalysisTaskSECompareHF
+// AliAnalysisTaskSE for the comparison of heavy-flavour decay candidates
+// to MC truth (kinematics stored in the AOD)
+// Author: A.Dainese, andrea.dainese@lnl.infn.it
+//*************************************************************************
+
+#include <TNtuple.h>
+#include <TH1F.h>
+
+#include "AliAnalysisTaskSE.h"
+
+class AliAnalysisTaskSECompareHF : public AliAnalysisTaskSE
+{
+ public:
+
+ AliAnalysisTaskSECompareHF();
+ AliAnalysisTaskSECompareHF(const char *name);
+ virtual ~AliAnalysisTaskSECompareHF();
+
+
+ // Implementation of interface methods
+ virtual void UserCreateOutputObjects();
+ virtual void Init();
+ virtual void LocalInit() {Init();}
+ virtual void UserExec(Option_t *option);
+ virtual void Terminate(Option_t *option);
+
+ void SetD0toKpiCuts(Double_t cut0=1000.,Double_t cut1=100000.,
+ Double_t cut2=1.1,Double_t cut3=0.,Double_t cut4=0.,
+ Double_t cut5=100000.,Double_t cut6=100000.,
+ Double_t cut7=100000000.,Double_t cut8=-1.1);
+ // cuts[0] = inv. mass half width [GeV]
+ // cuts[1] = dca [cm]
+ // cuts[2] = cosThetaStar
+ // cuts[3] = pTK [GeV/c]
+ // cuts[4] = pTPi [GeV/c]
+ // cuts[5] = d0K [cm] upper limit!
+ // cuts[6] = d0Pi [cm] upper limit!
+ // cuts[7] = d0d0 [cm^2]
+ // cuts[8] = cosThetaPoint
+ void SetD0toKpiCuts(const Double_t cuts[9]);
+
+ private:
+
+ AliAnalysisTaskSECompareHF(const AliAnalysisTaskSECompareHF &source);
+ AliAnalysisTaskSECompareHF& operator=(const AliAnalysisTaskSECompareHF& source);
+ TList *fOutput; //! list send on output slot 0
+ TNtuple *fNtupleD0Cmp; // output ntuple
+ TH1F *fHistMass; // output histogram
+ Double_t fD0toKpiCuts[9]; // cuts for D0->Kpi selection
+
+ ClassDef(AliAnalysisTaskSECompareHF,1); // AliAnalysisTaskSE for the MC association of heavy-flavour decay candidates
+};
+
+#endif
+
--- /dev/null
+void AliAnalysisTaskSECompareHFTest()\r
+{\r
+ //\r
+ // Test macro for the AliAnalysisTaskSE for heavy-flavour candidates\r
+ // association with MC truth (using MC info in AOD)\r
+ // A.Dainese, andrea.dainese@lnl.infn.it\r
+ //\r
+\r
+ gSystem->Load("libTree.so");\r
+ gSystem->Load("libGeom.so");\r
+ gSystem->Load("libPhysics.so");\r
+ gSystem->Load("libVMC.so");\r
+ gSystem->Load("libSTEERBase.so");\r
+ gSystem->Load("libESD.so");\r
+ gSystem->Load("libAOD.so"); \r
+ gSystem->Load("libANALYSIS.so");\r
+ gSystem->Load("libANALYSISalice.so");\r
+ gSystem->Load("libPWG3base.so");\r
+ gSystem->Load("libPWG3vertexingHF.so");\r
+\r
+\r
+ // Local files \r
+ TChain *chain = new TChain("aodTree");\r
+ chain->Add("./AliAOD.root");\r
+\r
+ // or:\r
+ /*\r
+ //Fetch files with AliEn :\r
+ const char *collectionfile = "CollectionTags.xml";\r
+ TGrid::Connect("alien://") ;\r
+ //Create an AliRunTagCuts and an AliEventTagCuts Object and impose some selection criteria\r
+ AliRunTagCuts *runCuts = new AliRunTagCuts();\r
+ AliEventTagCuts *eventCuts = new AliEventTagCuts();\r
+ AliLHCTagCuts *lhcCuts = new AliLHCTagCuts();\r
+ AliDetectorTagCuts *detCuts = new AliDetectorTagCuts();\r
+ eventCuts->SetMultiplicityRange(0,20000);\r
+ //Create an AliTagAnalysis Object and chain the tags\r
+ AliTagAnalysis *tagAna = new AliTagAnalysis();\r
+ tagAna->SetType("AOD");\r
+ TAlienCollection *coll = TAlienCollection::Open(collectionfile);\r
+ TGridResult *tagResult = coll->GetGridResult("",0,0);\r
+ tagResult->Print();\r
+ tagAna->ChainGridTags(tagResult);\r
+ //Create a new aod chain and assign the chain that is returned by querying the tags\r
+ TChain* chain = tagAna->QueryTags(runCuts,lhcCuts,detCuts,eventCuts);\r
+ */\r
+\r
+\r
+ // Create the analysis manager\r
+ AliAnalysisManager *mgr = new AliAnalysisManager("My Manager","My Manager");\r
+ mgr->SetDebugLevel(10);\r
+ \r
+\r
+ // Input\r
+ AliAODInputHandler *inputHandler = new AliAODInputHandler();\r
+ inputHandler->AddFriend("AliAOD.VertexingHF.root");\r
+ mgr->SetInputEventHandler(inputHandler);\r
+\r
+ \r
+ // Aanalysis task \r
+ AliAnalysisTaskSECompareHF *hfTask = new AliAnalysisTaskSECompareHF("CompareHFAnalysis");\r
+ hfTask->SetDebugLevel(2);\r
+ Double_t cutsD0[9]=\r
+ // cutsD0[0] = inv. mass half width [GeV]\r
+ // cutsD0[1] = dca [cm]\r
+ // cutsD0[2] = cosThetaStar\r
+ // cutsD0[3] = pTK [GeV/c]\r
+ // cutsD0[4] = pTPi [GeV/c]\r
+ // cutsD0[5] = d0K [cm] upper limit!\r
+ // cutsD0[6] = d0Pi [cm] upper limit!\r
+ // cutsD0[7] = d0d0 [cm^2]\r
+ // cutsD0[8] = cosThetaPoint\r
+ {1,\r
+ 100000.,\r
+ 1.1,\r
+ 0.,\r
+ 0.,\r
+ 100000.,\r
+ 100000.,\r
+ 100000000.,\r
+ -0.9}; \r
+ hfTask->SetD0toKpiCuts(cutsD0);\r
+ \r
+ mgr->AddTask(hfTask);\r
+ \r
+ //\r
+ // Create containers for input/output\r
+ AliAnalysisDataContainer *cinput = mgr->CreateContainer("cinput",TChain::Class(), \r
+ AliAnalysisManager::kInputContainer);\r
+ AliAnalysisDataContainer *coutput = mgr->CreateContainer("coutput",TList::Class(),\r
+ AliAnalysisManager::kOutputContainer, \r
+ "CmpHF.root");\r
+ mgr->ConnectInput(hfTask,0,cinput);\r
+ mgr->ConnectOutput(hfTask,1,coutput);\r
+\r
+ //\r
+ // Run the analysis\r
+ // \r
+ printf("CHAIN HAS %d ENTRIES\n",(Int_t)chain->GetEntries());\r
+ if(mgr->InitAnalysis()) {\r
+ mgr->PrintStatus();\r
+ mgr->StartAnalysis("local",chain);\r
+ //mgr->StartAnalysis("grid",chain);\r
+ }\r
+\r
+ return;\r
+}\r
git://git.uio.no / u/mrichter/AliRoot.git / commitdiff