PWG3/vertexingHF/AliAnalysisTaskSEDplus.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 extraction of signal(e.g D+) of heavy flavor
  19 // decay candidates with the MC truth.
  20 // Author: Renu Bala, bala@to.infn.it
  21 /////////////////////////////////////////////////////////////
  22
  23 #include <TClonesArray.h>
  24 #include <TNtuple.h>
  25 #include <TList.h>
  26 #include <TH1F.h>
  27 #include "AliAODEvent.h"
  28 #include "AliAODVertex.h"
  29 #include "AliAODTrack.h"
  30 #include "AliAODMCHeader.h"
  31 #include "AliAODMCParticle.h"
  32 #include "AliAODRecoDecayHF3Prong.h"
  33 #include "AliAnalysisVertexingHF.h"
  34 #include "AliAnalysisTaskSE.h"
  35 #include "AliAnalysisTaskSEDplus.h"
  36
  37 ClassImp(AliAnalysisTaskSEDplus)
  38
  39
  40 //________________________________________________________________________
  41 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
  42 AliAnalysisTaskSE(),
  43 fOutput(0),
  44 fNtupleDplus(0),
  45 fNtupleDplusbackg(0),
  46 fHistMass(0),
  47 fHistSignal(0),
  48 fHistBackground(0),
  49 fVHF(0)
  50 {
  51   // Default constructor
  52 }
  53
  54 //________________________________________________________________________
  55 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name):
  56 AliAnalysisTaskSE(name),
  57 fOutput(0),
  58 fNtupleDplus(0),
  59 fNtupleDplusbackg(0),
  60 fHistMass(0),
  61 fHistSignal(0),
  62 fHistBackground(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   // Output slot #2 writes into a TNtuple container
  70   DefineOutput(2,TNtuple::Class());  //My private output
  71   // Output slot #3 writes into a TNtuple container
  72   DefineOutput(3,TNtuple::Class());  //My private output
  73 }
  74
  75 //________________________________________________________________________
  76 AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
  77 {
  78   // Destructor
  79   if (fOutput) {
  80     delete fOutput;
  81     fOutput = 0;
  82   }
  83   if (fVHF) {
  84     delete fVHF;
  85     fVHF = 0;
  86   }
  87 }
  88
  89 //________________________________________________________________________
  90 void AliAnalysisTaskSEDplus::Init()
  91 {
  92   // Initialization
  93
  94   if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
  95
  96   gROOT->LoadMacro("ConfigVertexingHF.C");
  97
  98   fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
  99   fVHF->PrintStatus();
 100
 101   return;
 102 }
 103
 104 //________________________________________________________________________
 105 void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
 106 {
 107   // Create the output container
 108   //
 109   if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");
 110
 111   // Several histograms are more conveniently managed in a TList
 112   fOutput = new TList();
 113   fOutput->SetOwner();
 114
 115   fHistMass = new TH1F("fHistMass", "D^{+} invariant mass; M [GeV]; Entries",100,1.765,1.965);
 116   fHistSignal = new TH1F("fHistSignal", "D^{+} invariant mass - MC; M [GeV]; Entries",100,1.765,1.965);
 117   fHistBackground =new TH1F("fHistBackground", "Background invariant mass - MC; M [GeV]; Entries",100,1.765,1.965);
 118
 119
 120
 121
 122   fHistMass->Sumw2();
 123   fHistSignal->Sumw2();
 124   fHistBackground->Sumw2();
 125
 126
 127   //fHistMass->SetMinimum(0);
 128   //fHistSignal->SetMinimum(0);
 129   //fHistBackground->SetMinimum(0);
 130   fOutput->Add(fHistMass);
 131   fOutput->Add(fHistSignal);
 132   fOutput->Add(fHistBackground);
 133
 134   OpenFile(2); // 2 is the slot number of the ntuple
 135   fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:Ptpi:PtK:Ptpi2:PtRec:PtTrue:PointingAngle:DecLeng:VxTrue:VxRec:VyRec,VzRec,InvMass:sigvert:d0Pi:d0K:d0Pi2");
 136   OpenFile(3); // 3 is the slot number of the ntuple
 137   fNtupleDplusbackg = new TNtuple("fNtupleDplusbackg","D + backg","Ptpibkg:Ptkbkg:Ptpi2bkg:PtRecbkg:PointingAnglebkg:DLbkg:VxRecbkg:VyRecbkg:VzRecbkg:InvMassbkg:sigvertbkg:d0Pibkg:d0Kbkg:d0Pi2bkg");
 138
 139   return;
 140 }
 141
 142 //________________________________________________________________________
 143 void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
 144 {
 145   // Execute analysis for current event:
 146   // heavy flavor candidates association to MC truth
 147
 148
 149
 150   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 151
 152   // load Dplus->Kpipi candidates
 153   TClonesArray *array3Prong =
 154     (TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 155   if(!array3Prong) {
 156     printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
 157     return;
 158   }
 159
 160   // AOD primary vertex
 161     AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 162     //    vtx1->Print();
 163
 164   // load MC particles
 165   TClonesArray *arrayMC =
 166     (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 167   if(!arrayMC) {
 168     printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
 169     return;
 170   }
 171
 172   // load MC header
 173   AliAODMCHeader *mcHeader =
 174     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 175   if(!mcHeader) {
 176     printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
 177     return;
 178   }
 179     Int_t n3Prong = array3Prong->GetEntriesFast();
 180     if(fDebug>1) printf("Number of D+->Kpipi: %d\n",n3Prong);
 181
 182
 183     for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 184       AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 185
 186
 187       Bool_t unsetvtx=kFALSE;
 188       if(!d->GetOwnPrimaryVtx()){
 189         d->SetOwnPrimaryVtx(vtx1);
 190         unsetvtx=kTRUE;
 191       }
 192       if(d->SelectDplus(fVHF->GetDplusCuts())) {
 193
 194
 195         Int_t labDp = d->MatchToMC(411,arrayMC);
 196
 197         if(labDp>=0) {
 198           AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
 199           Int_t pdgDp = TMath::Abs(partDp->GetPdgCode());
 200           if(pdgDp==411){
 201
 202             fHistSignal->Fill(d->InvMassDplus());
 203             fHistMass->Fill(d->InvMassDplus());
 204
 205             // Post the data already here
 206             PostData(1,fOutput);
 207
 208             AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
 209
 210             Float_t tmp[20];
 211             tmp[0]=pdgDp;
 212             tmp[1]=partDp->Px()-d->Px();
 213             tmp[2]=partDp->Py()-d->Py();
 214             tmp[3]=partDp->Pz()-d->Pz();
 215             tmp[4]=d->PtProng(0);
 216             tmp[5]=d->PtProng(1);
 217             tmp[6]=d->PtProng(2);
 218             tmp[7]=d->Pt();
 219             tmp[8]=partDp->Pt();
 220             tmp[9]=d->CosPointingAngle();
 221             tmp[10]=d->DecayLength();
 222             tmp[11]=dg0->Xv();
 223             tmp[12]=d->Xv();
 224             tmp[13]=d->Yv();
 225             tmp[14]=d->Zv();
 226             tmp[15]=d->InvMassDplus();
 227             tmp[16]=d->GetSigmaVert();
 228             tmp[17]=d->Getd0Prong(0);
 229             tmp[18]=d->Getd0Prong(1);
 230             tmp[19]=d->Getd0Prong(2);
 231
 232
 233             fNtupleDplus->Fill(tmp);
 234             PostData(2,fNtupleDplus);
 235           }
 236         } else {
 237
 238             Float_t tmpbkg[14];
 239             tmpbkg[0]=d->PtProng(0);
 240             tmpbkg[1]=d->PtProng(1);
 241             tmpbkg[2]=d->PtProng(2);
 242             tmpbkg[3]=d->Pt();
 243             tmpbkg[4]=d->CosPointingAngle();
 244             tmpbkg[5]=d->DecayLength();
 245             tmpbkg[6]=d->Xv();
 246             tmpbkg[7]=d->Yv();
 247             tmpbkg[8]=d->Zv();
 248             tmpbkg[9]=d->InvMassDplus();
 249             tmpbkg[10]=d->GetSigmaVert();
 250             tmpbkg[11]=d->Getd0Prong(0);
 251             tmpbkg[12]=d->Getd0Prong(1);
 252             tmpbkg[13]=d->Getd0Prong(2);
 253
 254           fHistBackground->Fill(d->InvMassDplus());
 255           fNtupleDplusbackg->Fill(tmpbkg);
 256           PostData(3,fNtupleDplusbackg);
 257           fHistMass->Fill(d->InvMassDplus());
 258         }
 259
 260
 261       }
 262
 263       if(unsetvtx) d->UnsetOwnPrimaryVtx();
 264
 265     }
 266      // end loop on D+->Kpipi
 267
 268
 269   return;
 270 }
 271
 272 //________________________________________________________________________
 273 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
 274 {
 275   // Terminate analysis
 276   //
 277   if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
 278
 279   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 280   if (!fOutput) {
 281     printf("ERROR: fOutput not available\n");
 282     return;
 283   }
 284
 285   fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
 286   fHistSignal = dynamic_cast<TH1F*>(fOutput->FindObject("fHistSignal"));
 287   fHistBackground = dynamic_cast<TH1F*>(fOutput->FindObject("fHistBackground"));
 288   fNtupleDplus = dynamic_cast<TNtuple*>(GetOutputData(2));
 289   fNtupleDplusbackg = dynamic_cast<TNtuple*>(GetOutputData(3));
 290
 291   return;
 292 }
 293