PWG3/vertexingHF/AliAnalysisTaskSECompareHF.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /////////////////////////////////////////////////////////////
  17 //
  18 // AliAnalysisTaskSE for the comparison of heavy flavor
  19 // decay candidates with the MC truth.
  20 //
  21 // Author: A.Dainese, andrea.dainese@lnl.infn.it
  22 /////////////////////////////////////////////////////////////
  23
  24 #include <TClonesArray.h>
  25 #include <TNtuple.h>
  26 #include <TList.h>
  27 #include <TH1F.h>
  28
  29 #include "AliAODEvent.h"
  30 #include "AliAODVertex.h"
  31 #include "AliAODTrack.h"
  32 #include "AliAODMCHeader.h"
  33 #include "AliAODMCParticle.h"
  34 #include "AliAODRecoDecayHF2Prong.h"
  35 #include "AliAODRecoCascadeHF.h"
  36 #include "AliAnalysisVertexingHF.h"
  37 #include "AliAnalysisTaskSE.h"
  38 #include "AliAnalysisTaskSECompareHF.h"
  39
  40 ClassImp(AliAnalysisTaskSECompareHF)
  41
  42
  43 //________________________________________________________________________
  44 AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF():
  45 AliAnalysisTaskSE(),
  46 fOutput(0),
  47 fNtupleD0Cmp(0),
  48 fHistMass(0),
  49 fVHF(0)
  50 {
  51   // Default constructor
  52
  53   // Output slot #1 writes into a TList container
  54   DefineOutput(1,TList::Class());  //My private output
  55 }
  56
  57 //________________________________________________________________________
  58 AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name):
  59 AliAnalysisTaskSE(name),
  60 fOutput(0),
  61 fNtupleD0Cmp(0),
  62 fHistMass(0),
  63 fVHF(0)
  64 {
  65   // Default constructor
  66
  67   // Output slot #1 writes into a TList container
  68   DefineOutput(1,TList::Class());  //My private output
  69 }
  70
  71 //________________________________________________________________________
  72 AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF()
  73 {
  74   // Destructor
  75   if (fOutput) {
  76     delete fOutput;
  77     fOutput = 0;
  78   }
  79   if (fVHF) {
  80     delete fVHF;
  81     fVHF = 0;
  82   }
  83 }
  84
  85 //________________________________________________________________________
  86 void AliAnalysisTaskSECompareHF::Init()
  87 {
  88   // Initialization
  89
  90   if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n");
  91
  92   gROOT->LoadMacro("ConfigVertexingHF.C");
  93
  94   fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
  95   fVHF->PrintStatus();
  96
  97   return;
  98 }
  99
 100 //________________________________________________________________________
 101 void AliAnalysisTaskSECompareHF::UserCreateOutputObjects()
 102 {
 103   // Create the output container
 104   //
 105   if(fDebug > 1) printf("AnalysisTaskSECompareHF::UserCreateOutputObjects() \n");
 106
 107   // Several histograms are more conveniently managed in a TList
 108   fOutput = new TList();
 109   fOutput->SetOwner();
 110
 111   fHistMass = new TH1F("fHistMass", "D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
 112   fHistMass->Sumw2();
 113   fHistMass->SetMinimum(0);
 114   fOutput->Add(fHistMass);
 115
 116   fNtupleD0Cmp = new TNtuple("fNtupleD0Cmp","D0 comparison","pdg:VxRec:VxTrue:PtRec:PtTrue");
 117   fOutput->Add(fNtupleD0Cmp);
 118
 119   return;
 120 }
 121
 122 //________________________________________________________________________
 123 void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
 124 {
 125   // Execute analysis for current event:
 126   // heavy flavor candidates association to MC truth
 127
 128   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 129
 130   // load D0->Kpi candidates
 131   TClonesArray *inputArrayD0toKpi =
 132     (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
 133   if(!inputArrayD0toKpi) {
 134     printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
 135     return;
 136   }
 137
 138   /*
 139   // load D*+ candidates
 140   TClonesArray *inputArrayDstar =
 141     (TClonesArray*)aod->GetList()->FindObject("Dstar");
 142   if(!inputArrayDstar) {
 143     printf("AliAnalysisTaskSECompareHF::UserExec: Dstar branch not found!\n");
 144     return;
 145   }
 146   */
 147
 148   // AOD primary vertex
 149   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 150   //vtx1->Print();
 151
 152   // load MC particles
 153   TClonesArray *mcArray =
 154     (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 155   if(!mcArray) {
 156     printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n");
 157     return;
 158   }
 159
 160   // load MC header
 161   AliAODMCHeader *mcHeader =
 162     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 163   if(!mcHeader) {
 164     printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n");
 165     return;
 166   }
 167
 168
 169   // loop over D0->Kpi candidates
 170   Int_t nInD0toKpi = inputArrayD0toKpi->GetEntriesFast();
 171   printf("Number of D0->Kpi: %d\n",nInD0toKpi);
 172
 173   for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
 174     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayD0toKpi->UncheckedAt(iD0toKpi);
 175     Bool_t unsetvtx=kFALSE;
 176     if(!d->GetOwnPrimaryVtx()) {
 177       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 178       unsetvtx=kTRUE;
 179     }
 180     Int_t okD0=0,okD0bar=0;
 181     if(d->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) {
 182       Int_t labD0 = d->MatchToMC(421,mcArray);
 183       if(labD0>=0) {
 184         AliAODMCParticle *partD0 = (AliAODMCParticle*)mcArray->At(labD0);
 185         Int_t pdgD0 = partD0->GetPdgCode();
 186         Double_t invmass = (pdgD0==421 ? d->InvMassD0() : d->InvMassD0bar());
 187         fHistMass->Fill(invmass);
 188         // Post the data already here
 189         PostData(1,fOutput);
 190
 191         fNtupleD0Cmp->Fill(pdgD0,d->Xv(),partD0->Xv(),d->Pt(),partD0->Pt());
 192       }
 193     }
 194     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 195   } // end loop on D0->Kpi
 196
 197
 198   // SIMPLE EXAMPLE OF D*+ MATCHING
 199   /*
 200   // loop over D*+ candidates
 201   for (Int_t iDstar = 0; iDstar < inputArrayDstar->GetEntries(); iDstar++) {
 202     AliAODRecoCascadeHF *c = (AliAODRecoCascadeHF*)inputArrayDstar->UncheckedAt(iDstar);
 203     Int_t labDstar = c->MatchToMC(413,421,mcArray);
 204     if(labDstar>=0) printf("GOOD MATCH FOR D*+\n");
 205   }
 206   */
 207
 208   return;
 209 }
 210
 211 //________________________________________________________________________
 212 void AliAnalysisTaskSECompareHF::Terminate(Option_t */*option*/)
 213 {
 214   // Terminate analysis
 215   //
 216   if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n");
 217
 218   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 219   if (!fOutput) {
 220     printf("ERROR: fOutput not available\n");
 221     return;
 222   }
 223
 224   fNtupleD0Cmp = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleD0Cmp"));
 225   fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
 226
 227   return;
 228 }
 229