PWG3/vertexingHF/AliAnalysisTaskSED0Mass.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, 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 D0 candidates invariant mass histogram
  19 // and comparison with the MC truth.
  20 //
  21 // Authors: A.Dainese, andrea.dainese@lnl.infn.it
  22 // and Chiara Bianchin, chiara.bianchin@pd.infn.it
  23 /////////////////////////////////////////////////////////////
  24
  25 #include <Riostream.h>
  26 #include <TClonesArray.h>
  27 #include <TNtuple.h>
  28 #include <TList.h>
  29 #include <TH1F.h>
  30
  31
  32 #include "AliAODEvent.h"
  33 #include "AliAODVertex.h"
  34 #include "AliAODTrack.h"
  35 #include "AliAODMCHeader.h"
  36 #include "AliAODMCParticle.h"
  37 #include "AliAODRecoDecayHF2Prong.h"
  38 #include "AliAODRecoCascadeHF.h"
  39 #include "AliAnalysisVertexingHF.h"
  40 #include "AliAnalysisTaskSE.h"
  41 #include "AliAnalysisTaskSED0Mass.h"
  42
  43 ClassImp(AliAnalysisTaskSED0Mass)
  44
  45
  46 //________________________________________________________________________
  47 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass():
  48 AliAnalysisTaskSE(),
  49 fOutput(0),
  50 fVHF(0)
  51 {
  52   // Default constructor
  53 }
  54
  55 //________________________________________________________________________
  56 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name):
  57 AliAnalysisTaskSE(name),
  58 fOutput(0),
  59 fVHF(0)
  60 {
  61   // Default constructor
  62
  63   // Output slot #1 writes into a TList container
  64   DefineOutput(1,TList::Class());  //My private output
  65 }
  66
  67 //________________________________________________________________________
  68 AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass()
  69 {
  70   // Destructor
  71   if (fOutput) {
  72     delete fOutput;
  73     fOutput = 0;
  74   }
  75   if (fVHF) {
  76     delete fVHF;
  77     fVHF = 0;
  78   }
  79
  80 }
  81
  82 //________________________________________________________________________
  83 void AliAnalysisTaskSED0Mass::Init()
  84 {
  85   // Initialization
  86
  87   if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
  88
  89   gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/ConfigVertexingHF.C");
  90
  91   fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
  92   // set dedidcated cuts
  93   fVHF->SetD0toKpiCuts(0.7,0.03,0.8,0.06,0.06,0.05,0.05,-0.0002,0.6); //2.p-p vertex reconstructed
  94   fVHF->PrintStatus();
  95
  96   return;
  97 }
  98
  99 //________________________________________________________________________
 100 void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
 101 {
 102
 103   // Create the output container
 104   //
 105   if(fDebug > 1) printf("AnalysisTaskSED0Mass::UserCreateOutputObjects() \n");
 106
 107   // Several histograms are more conveniently managed in a TList
 108   fOutput = new TList();
 109   fOutput->SetOwner();
 110
 111   const Int_t nhist=3;
 112   TString nameMass, nameSgn, nameBkg;
 113
 114   for(Int_t i=0;i<nhist;i++){
 115     nameMass="histMass_";
 116     nameMass+=i+1;
 117     TH1F* tmpM = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
 118     nameSgn="histSgn_";
 119     nameSgn+=i+1;
 120     TH1F* tmpS = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
 121     nameBkg="histBkg_";
 122     nameBkg+=i+1;
 123     TH1F* tmpB = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
 124     printf("Created histograms %s\t%s\t%s\n",tmpM->GetName(),tmpS->GetName(),tmpB->GetName());
 125
 126     fOutput->Add(tmpM);
 127     fOutput->Add(tmpS);
 128     fOutput->Add(tmpB);
 129
 130   }
 131   fOutput->ls();
 132   return;
 133 }
 134
 135 //________________________________________________________________________
 136 void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
 137 {
 138   // Execute analysis for current event:
 139   // heavy flavor candidates association to MC truth
 140
 141   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 142
 143   // load D0->Kpi candidates
 144   TClonesArray *inputArrayD0toKpi =
 145     (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
 146   if(!inputArrayD0toKpi) {
 147     printf("AliAnalysisTaskSECompareHFpt::UserExec: D0toKpi branch not found!\n");
 148     return;
 149   }
 150
 151   // AOD primary vertex
 152   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 153   //vtx1->Print();
 154
 155   // load MC particles
 156   TClonesArray *mcArray =
 157     (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 158   if(!mcArray) {
 159     printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
 160     return;
 161   }
 162
 163   // load MC header
 164   AliAODMCHeader *mcHeader =
 165     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 166   if(!mcHeader) {
 167     printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
 168     return;
 169   }
 170
 171
 172   // loop over D0->Kpi candidates
 173   Int_t nInD0toKpi = inputArrayD0toKpi->GetEntriesFast();
 174   printf("Number of D0->Kpi: %d\n",nInD0toKpi);
 175
 176   for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
 177
 178     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayD0toKpi->UncheckedAt(iD0toKpi);
 179     Bool_t unsetvtx=kFALSE;
 180     if(!d->GetOwnPrimaryVtx()) {
 181       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 182       unsetvtx=kTRUE;
 183     }
 184
 185     //  cout<<"Cuts applied: ";
 186     //     Double_t *cutsapplied=(Double_t*)fVHF->GetD0toKpiCuts();
 187     //     for (Int_t i=0;i<9;i++){
 188     //       printf(cutsapplied[i]\t);
 189     //     }
 190     //    cout<<endl;
 191
 192
 193     Double_t pt = d->Pt();
 194     Int_t ptbin=0;
 195
 196     //cout<<"P_t = "<<pt<<endl;
 197     if (pt>0. && pt<=1.) {
 198       ptbin=1;
 199       //printf("I'm in the bin %d\n",ptbin);
 200     }
 201     else {
 202       if(pt>1. && pt<=3.) {
 203         ptbin=2;
 204         //printf("I'm in the bin %d\n",ptbin);
 205       }
 206       else {
 207         ptbin=3;
 208         //printf("I'm in the bin %d\n",ptbin);
 209       }//if(pt>3)
 210     }
 211
 212     FillHists(ptbin,d,mcArray);
 213
 214     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 215
 216   }
 217
 218
 219
 220   // Post the data
 221   PostData(1,fOutput);
 222
 223
 224   return;
 225 }
 226 //____________________________________________________________________________*
 227 void AliAnalysisTaskSED0Mass::FillHists(Int_t ptbin, AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC){
 228   //
 229   // function used in UserExec:
 230   //
 231   Int_t okD0=0,okD0bar=0;
 232
 233   if(part->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) {//selected
 234     Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
 235     //printf("SELECTED\n");
 236     Int_t labD0 = part->MatchToMC(421,arrMC); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
 237     //printf("labD0 %d",labD0);
 238
 239     if(labD0>=0) {
 240
 241       AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
 242
 243       Int_t pdgD0 = partD0->GetPdgCode();
 244
 245       //printf(" pdgD0 %d\n",pdgD0);
 246       TString fillthis="histSgn_";
 247       fillthis+=ptbin;
 248       cout<<"Filling "<<fillthis<<endl;
 249
 250       if (pdgD0==421){ //D0
 251         //cout<<"Fill S with D0"<<endl;
 252         ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(invmassD0);
 253         //cout<<"Address "<<fOutput->FindObject(fillthis)<<endl;
 254         //cout<<"Name "<<(fOutput->FindObject(fillthis))->GetName()<<endl;
 255       }
 256       else {//D0bar
 257         //printf("Fill S with D0bar");
 258         ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(invmassD0bar);
 259       }
 260
 261     }
 262     else {//background
 263       TString fillthis="histBkg_";
 264       fillthis+=ptbin;
 265       //cout<<"Filling "<<fillthis<<endl;
 266
 267           //printf("Fill background");
 268       ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(invmassD0);
 269       ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(invmassD0bar);
 270     }
 271
 272     //no MC info, just cut selection
 273     TString fillthis="histMass_";
 274     fillthis+=ptbin;
 275     cout<<"Filling "<<fillthis<<endl;
 276
 277     if (okD0==1) {
 278       //printf("Fill mass with D0");
 279       ((TH1F*)(fOutput->FindObject(fillthis)))->Fill(invmassD0);
 280
 281     }
 282     if (okD0bar==1) {
 283       ((TH1F*)fOutput->FindObject(fillthis))->Fill(invmassD0bar);
 284       //printf("Fill mass with D0bar");
 285     }
 286
 287   } //else cout<<"NOT SELECTED"<<endl;
 288
 289 }
 290 //________________________________________________________________________
 291 void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
 292 {
 293   // Terminate analysis
 294   //
 295   if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
 296
 297   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 298   if (!fOutput) {
 299     printf("ERROR: fOutput not available\n");
 300     return;
 301   }
 302
 303   return;
 304 }
 305