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 <TString.h>
  27 #include <TH1F.h>
  28 #include "AliAODEvent.h"
  29 #include "AliAODVertex.h"
  30 #include "AliAODTrack.h"
  31 #include "AliAODMCHeader.h"
  32 #include "AliAODMCParticle.h"
  33 #include "AliAODRecoDecayHF3Prong.h"
  34 #include "AliAnalysisVertexingHF.h"
  35 #include "AliAnalysisTaskSE.h"
  36 #include "AliAnalysisTaskSEDplus.h"
  37
  38 ClassImp(AliAnalysisTaskSEDplus)
  39
  40
  41 //________________________________________________________________________
  42 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
  43 AliAnalysisTaskSE(),
  44 fOutput(0),
  45 fNtupleDplus(0),
  46 fNtupleDplusbackg(0),
  47 fFillNtuple(kFALSE),
  48 fVHF(0)
  49 {
  50   // Default constructor
  51 }
  52
  53 //________________________________________________________________________
  54 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,Bool_t fillNtuple):
  55 AliAnalysisTaskSE(name),
  56 fOutput(0),
  57 fNtupleDplus(0),
  58 fNtupleDplusbackg(0),
  59 fFillNtuple(fillNtuple),
  60 fVHF(0)
  61 {
  62   // Default constructor
  63
  64   // Output slot #1 writes into a TList container
  65   DefineOutput(1,TList::Class());  //My private output
  66
  67   if(fFillNtuple){
  68     // Output slot #2 writes into a TNtuple container
  69     DefineOutput(2,TNtuple::Class());  //My private output
  70     // Output slot #3 writes into a TNtuple container
  71     DefineOutput(3,TNtuple::Class());  //My private output
  72   }
  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   fOutput->SetName("OutputHistos");
 115
 116   Int_t nPtBins=4;
 117   TString hisname;
 118   for(Int_t i=0;i<nPtBins;i++){
 119     hisname.Form("hMassPt%d",i);
 120     TH1F* hm=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
 121     hm->Sumw2();
 122     hisname.Form("hSigPt%d",i);
 123     TH1F* hs=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
 124     hs->Sumw2();
 125     hisname.Form("hBkgPt%d",i);
 126     TH1F* hb=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
 127     hb->Sumw2();
 128
 129     hisname.Form("hMassPt%dTC",i);
 130     TH1F* hmtc=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
 131     hmtc->Sumw2();
 132     hisname.Form("hSigPt%dTC",i);
 133     TH1F* hstc=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
 134     hstc->Sumw2();
 135     hisname.Form("hBkgPt%dTC",i);
 136     TH1F* hbtc=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
 137     hbtc->Sumw2();
 138
 139
 140     fOutput->Add(hm);
 141     fOutput->Add(hs);
 142     fOutput->Add(hb);
 143     fOutput->Add(hmtc);
 144     fOutput->Add(hstc);
 145     fOutput->Add(hbtc);
 146   }
 147
 148
 149   if(fFillNtuple){
 150     OpenFile(2); // 2 is the slot number of the ntuple
 151     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");
 152     OpenFile(3); // 3 is the slot number of the ntuple
 153     fNtupleDplusbackg = new TNtuple("fNtupleDplusbackg","D + backg","Ptpibkg:Ptkbkg:Ptpi2bkg:PtRecbkg:PointingAnglebkg:DLbkg:VxRecbkg:VyRecbkg:VzRecbkg:InvMassbkg:sigvertbkg:d0Pibkg:d0Kbkg:d0Pi2bkg");
 154   }
 155
 156   return;
 157 }
 158
 159 //________________________________________________________________________
 160 void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
 161 {
 162   // Execute analysis for current event:
 163   // heavy flavor candidates association to MC truth
 164
 165
 166
 167   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 168
 169   // In case there is an AOD handler writing a standard AOD, use the AOD
 170   // event in memory rather than the input (ESD) event.
 171   if (!aod && AODEvent() && IsStandardAOD()) aod = dynamic_cast<AliAODEvent*> (AODEvent());
 172
 173   // load Dplus->Kpipi candidates
 174   TClonesArray *array3Prong =
 175     (TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 176   if(!array3Prong) {
 177     printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
 178     return;
 179   }
 180
 181   // AOD primary vertex
 182   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 183   //    vtx1->Print();
 184
 185   // load MC particles
 186   TClonesArray *arrayMC =
 187     (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 188   if(!arrayMC) {
 189     printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
 190     return;
 191   }
 192
 193   // load MC header
 194   AliAODMCHeader *mcHeader =
 195     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 196   if(!mcHeader) {
 197     printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
 198     return;
 199   }
 200   Int_t n3Prong = array3Prong->GetEntriesFast();
 201   if(fDebug>1) printf("Number of D+->Kpipi: %d\n",n3Prong);
 202
 203
 204   Int_t pdgDgDplustoKpipi[3]={321,211,211};
 205   Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};
 206   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 207     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 208
 209
 210     Bool_t unsetvtx=kFALSE;
 211     if(!d->GetOwnPrimaryVtx()){
 212       d->SetOwnPrimaryVtx(vtx1);
 213       unsetvtx=kTRUE;
 214     }
 215     if(d->SelectDplus(fVHF->GetDplusCuts())) {
 216       Int_t iPtBin=0;
 217       Double_t ptCand = d->Pt();
 218       if(ptCand<2.){
 219         iPtBin=0;
 220         cutsDplus[7]=0.08;
 221         cutsDplus[8]=0.5;
 222         cutsDplus[10]=0.979;
 223       }
 224       else if(ptCand>2. && ptCand<3){
 225         iPtBin=1;
 226         cutsDplus[7]=0.08;
 227         cutsDplus[8]=0.5;
 228         cutsDplus[9]=0.991;
 229       }else if(ptCand>3. && ptCand<5){
 230         iPtBin=2;
 231         cutsDplus[7]=0.1;
 232         cutsDplus[8]=0.5;
 233         cutsDplus[9]=0.9955;
 234       }else{
 235         iPtBin=3;
 236         cutsDplus[7]=0.1;
 237         cutsDplus[8]=0.5;
 238         cutsDplus[9]=0.997;
 239       }
 240       Bool_t passTightCuts=d->SelectDplus(cutsDplus);
 241       Int_t labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
 242       Double_t invMass=d->InvMassDplus();
 243
 244       TString hisNameA(Form("hMassPt%d",iPtBin));
 245       TString hisNameS(Form("hSigPt%d",iPtBin));
 246       TString hisNameB(Form("hBkgPt%d",iPtBin));
 247       TString hisNameATC(Form("hMassPt%dTC",iPtBin));
 248       TString hisNameSTC(Form("hSigPt%dTC",iPtBin));
 249       TString hisNameBTC(Form("hBkgPt%dTC",iPtBin));
 250
 251       ((TH1F*)(fOutput->FindObject(hisNameA)))->Fill(invMass);
 252       if(passTightCuts){
 253         ((TH1F*)(fOutput->FindObject(hisNameATC)))->Fill(invMass);
 254       }
 255
 256
 257       if(labDp>=0) {
 258         ((TH1F*)(fOutput->FindObject(hisNameS)))->Fill(invMass);
 259         if(passTightCuts){
 260           ((TH1F*)(fOutput->FindObject(hisNameSTC)))->Fill(invMass);
 261         }
 262         if(fFillNtuple){
 263           AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
 264           AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
 265           Float_t tmp[20];
 266           tmp[0]=TMath::Abs(partDp->GetPdgCode());
 267           tmp[1]=partDp->Px()-d->Px();
 268           tmp[2]=partDp->Py()-d->Py();
 269           tmp[3]=partDp->Pz()-d->Pz();
 270           tmp[4]=d->PtProng(0);
 271           tmp[5]=d->PtProng(1);
 272           tmp[6]=d->PtProng(2);
 273           tmp[7]=d->Pt();
 274           tmp[8]=partDp->Pt();
 275           tmp[9]=d->CosPointingAngle();
 276           tmp[10]=d->DecayLength();
 277           tmp[11]=dg0->Xv();
 278           tmp[12]=d->Xv();
 279           tmp[13]=d->Yv();
 280           tmp[14]=d->Zv();
 281           tmp[15]=d->InvMassDplus();
 282           tmp[16]=d->GetSigmaVert();
 283           tmp[17]=d->Getd0Prong(0);
 284           tmp[18]=d->Getd0Prong(1);
 285           tmp[19]=d->Getd0Prong(2);
 286           fNtupleDplus->Fill(tmp);
 287           PostData(2,fNtupleDplus);
 288         }
 289       }else{
 290         ((TH1F*)(fOutput->FindObject(hisNameB)))->Fill(invMass);
 291         if(passTightCuts){
 292           ((TH1F*)(fOutput->FindObject(hisNameBTC)))->Fill(invMass);
 293         }
 294         if(fFillNtuple){
 295           Float_t tmpbkg[14];
 296           tmpbkg[0]=d->PtProng(0);
 297           tmpbkg[1]=d->PtProng(1);
 298           tmpbkg[2]=d->PtProng(2);
 299           tmpbkg[3]=d->Pt();
 300           tmpbkg[4]=d->CosPointingAngle();
 301           tmpbkg[5]=d->DecayLength();
 302           tmpbkg[6]=d->Xv();
 303           tmpbkg[7]=d->Yv();
 304           tmpbkg[8]=d->Zv();
 305           tmpbkg[9]=d->InvMassDplus();
 306           tmpbkg[10]=d->GetSigmaVert();
 307           tmpbkg[11]=d->Getd0Prong(0);
 308           tmpbkg[12]=d->Getd0Prong(1);
 309           tmpbkg[13]=d->Getd0Prong(2);
 310           fNtupleDplusbackg->Fill(tmpbkg);
 311           PostData(3,fNtupleDplusbackg);
 312         }
 313       }
 314       PostData(1,fOutput);
 315     }
 316     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 317   }
 318
 319   return;
 320 }
 321
 322 //________________________________________________________________________
 323 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
 324 {
 325   // Terminate analysis
 326   //
 327   if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
 328
 329   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 330   if (!fOutput) {
 331     printf("ERROR: fOutput not available\n");
 332     return;
 333   }
 334
 335   if(fFillNtuple){
 336     fNtupleDplus = dynamic_cast<TNtuple*>(GetOutputData(2));
 337     fNtupleDplusbackg = dynamic_cast<TNtuple*>(GetOutputData(3));
 338   }
 339
 340   return;
 341 }
 342