New class for AOD<->MC association

diff --git a/PWG3/vertexingHF/AliAnalysisTaskSECompareHF.cxx b/PWG3/vertexingHF/AliAnalysisTaskSECompareHF.cxx
new file mode 100644 (file)
index 0000000..1b44549
--- /dev/null
+++ b/PWG3/vertexingHF/AliAnalysisTaskSECompareHF.cxx
@@ -0,0 +1,301 @@
+/**************************************************************************
+ * 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;
+}

diff --git a/PWG3/vertexingHF/AliAnalysisTaskSECompareHF.h b/PWG3/vertexingHF/AliAnalysisTaskSECompareHF.h
new file mode 100644 (file)
index 0000000..bc9a48a
--- /dev/null
+++ b/PWG3/vertexingHF/AliAnalysisTaskSECompareHF.h
@@ -0,0 +1,63 @@
+#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
+

diff --git a/PWG3/vertexingHF/AliAnalysisTaskSECompareHFTest.C b/PWG3/vertexingHF/AliAnalysisTaskSECompareHFTest.C
new file mode 100644 (file)
index 0000000..d2431e5
--- /dev/null
+++ b/PWG3/vertexingHF/AliAnalysisTaskSECompareHFTest.C
@@ -0,0 +1,107 @@
+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