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 "AliAnalysisVertexingHF.h"
  36 #include "AliAnalysisTaskSE.h"
  37 #include "AliAnalysisTaskSECompareHF.h"
  38
  39 ClassImp(AliAnalysisTaskSECompareHF)
  40
  41
  42 //________________________________________________________________________
  43 AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF():
  44 AliAnalysisTaskSE(),
  45 fOutput(0),
  46 fNtupleD0Cmp(0),
  47 fHistMass(0),
  48 fVHF(0)
  49 {
  50   // Default constructor
  51
  52   // Output slot #1 writes into a TList container
  53   DefineOutput(1,TList::Class());  //My private output
  54 }
  55
  56 //________________________________________________________________________
  57 AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name):
  58 AliAnalysisTaskSE(name),
  59 fOutput(0),
  60 fNtupleD0Cmp(0),
  61 fHistMass(0),
  62 fVHF(0)
  63 {
  64   // Default constructor
  65
  66   // Output slot #1 writes into a TList container
  67   DefineOutput(1,TList::Class());  //My private output
  68 }
  69
  70 //________________________________________________________________________
  71 AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF()
  72 {
  73   // Destructor
  74   if (fOutput) {
  75     delete fOutput;
  76     fOutput = 0;
  77   }
  78   if (fVHF) {
  79     delete fVHF;
  80     fVHF = 0;
  81   }
  82 }
  83
  84 //________________________________________________________________________
  85 void AliAnalysisTaskSECompareHF::Init()
  86 {
  87   // Initialization
  88
  89   if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n");
  90
  91   gROOT->LoadMacro("ConfigVertexingHF.C");
  92
  93   fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
  94   fVHF->PrintStatus();
  95
  96   return;
  97 }
  98
  99 //________________________________________________________________________
 100 void AliAnalysisTaskSECompareHF::UserCreateOutputObjects()
 101 {
 102   // Create the output container
 103   //
 104   if(fDebug > 1) printf("AnalysisTaskSECompareHF::UserCreateOutputObjects() \n");
 105
 106   // Several histograms are more conveniently managed in a TList
 107   fOutput = new TList();
 108   fOutput->SetOwner();
 109
 110   fHistMass = new TH1F("fHistMass", "D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
 111   fHistMass->Sumw2();
 112   fHistMass->SetMinimum(0);
 113   fOutput->Add(fHistMass);
 114
 115   fNtupleD0Cmp = new TNtuple("fNtupleD0Cmp","D0 comparison","pdg:VxRec:VxTrue:PtRec:PtTrue");
 116   fOutput->Add(fNtupleD0Cmp);
 117
 118   return;
 119 }
 120
 121 //________________________________________________________________________
 122 void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
 123 {
 124   // Execute analysis for current event:
 125   // heavy flavor candidates association to MC truth
 126
 127   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 128
 129   // load D0->Kpi candidates
 130   TClonesArray *inputArrayD0toKpi =
 131     (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
 132   if(!inputArrayD0toKpi) {
 133     printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
 134     return;
 135   }
 136
 137   // AOD primary vertex
 138   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 139   //vtx1->Print();
 140
 141   // load MC particles
 142   TClonesArray *mcArray =
 143     (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 144   if(!mcArray) {
 145     printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n");
 146     return;
 147   }
 148
 149   // load MC header
 150   AliAODMCHeader *mcHeader =
 151     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 152   if(!mcHeader) {
 153     printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n");
 154     return;
 155   }
 156
 157
 158
 159   // loop over D0->Kpi candidates
 160   Int_t nInD0toKpi = inputArrayD0toKpi->GetEntriesFast();
 161   printf("Number of D0->Kpi: %d\n",nInD0toKpi);
 162
 163   Int_t lab0,lab1,labMother,labD0daugh0,labD0daugh1,pdgMother,pdgD0;
 164
 165   for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
 166     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayD0toKpi->UncheckedAt(iD0toKpi);
 167     labD0daugh0=-1;
 168     labD0daugh1=-1;
 169     Bool_t unsetvtx=kFALSE;
 170     if(!d->GetOwnPrimaryVtx()) {
 171       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 172       unsetvtx=kTRUE;
 173     }
 174     Int_t okD0=0,okD0bar=0;
 175     if(d->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) {
 176       // get daughter AOD tracks
 177       AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
 178       AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1);
 179
 180       lab0 = trk0->GetLabel();
 181       lab1 = trk1->GetLabel();
 182
 183
 184       AliAODMCParticle *part0 = (AliAODMCParticle*)mcArray->At(lab0);
 185       if(!part0) {
 186         printf("no MC particle\n");
 187         continue;
 188       }
 189       while(part0->GetMother()>=0) {
 190         labMother=part0->GetMother();
 191         part0 = (AliAODMCParticle*)mcArray->At(labMother);
 192         if(!part0) {
 193           printf("no MC mother particle\n");
 194           break;
 195         }
 196         pdgMother = TMath::Abs(part0->GetPdgCode());
 197         if(pdgMother==421) {
 198           labD0daugh0=labMother;
 199           break;
 200         }
 201       }
 202
 203       AliAODMCParticle *part1 = (AliAODMCParticle*)mcArray->At(lab1);
 204       if(!part1) {
 205         printf("no MC particle\n");
 206         continue;
 207       }
 208       while(part1->GetMother()>=0) {
 209         labMother=part1->GetMother();
 210         part1 = (AliAODMCParticle*)mcArray->At(labMother);
 211         if(!part1) {
 212           printf("no MC mother particle\n");
 213           break;
 214         }
 215         pdgMother = TMath::Abs(part1->GetPdgCode());
 216         if(pdgMother==421) {
 217           labD0daugh1=labMother;
 218           break;
 219         }
 220       }
 221
 222       // check if the candidate is signal
 223       if(labD0daugh0>=0 && labD0daugh1>=0 && labD0daugh0==labD0daugh1) {
 224
 225         AliAODMCParticle *partD0 = (AliAODMCParticle*)mcArray->At(labD0daugh0);
 226         // check that the D0 decays in 2 prongs
 227         if (TMath::Abs(partD0->GetDaughter(1)-partD0->GetDaughter(0))==1) {
 228
 229           pdgD0 = partD0->GetPdgCode();
 230           Double_t invmass = (pdgD0==421 ? d->InvMassD0() : d->InvMassD0bar());
 231           fHistMass->Fill(invmass);
 232
 233           // Post the data already here
 234           PostData(1,fOutput);
 235
 236           fNtupleD0Cmp->Fill(pdgD0,d->Xv(),partD0->Xv(),d->Pt(),partD0->Pt());
 237         }
 238       }
 239
 240
 241     }
 242     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 243   } // end loop on D0->Kpi
 244
 245
 246   return;
 247 }
 248
 249 //________________________________________________________________________
 250 void AliAnalysisTaskSECompareHF::Terminate(Option_t */*option*/)
 251 {
 252   // Terminate analysis
 253   //
 254   if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n");
 255
 256   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 257   if (!fOutput) {
 258     printf("ERROR: fOutput not available\n");
 259     return;
 260   }
 261
 262   fNtupleD0Cmp = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleD0Cmp"));
 263   fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
 264
 265   return;
 266 }
 267