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