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