PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.cxx

   1 #include "AliAnalysisTaskSE.h"
   2 #include "AliAnalysisManager.h"
   3 #include "AliAnalysisDataContainer.h"
   4 #include "AliESDEvent.h"
   5 #include "AliESDtrackCuts.h"
   6 #include "AliStack.h"
   7 #include "AliCentrality.h"
   8 #include "AliMCEventHandler.h"
   9 #include "AliMCEvent.h"
  10 #include "AliHeader.h"
  11 #include "AliGenCocktailEventHeader.h"
  12 #include "AliGenEventHeader.h"
  13 #include "AliGenerator.h"
  14 #include "AliVertexingHFUtils.h"
  15 #include "AliMultiplicity.h"
  16 #include <TParticle.h>
  17 #include <TSystem.h>
  18 #include <TTree.h>
  19 #include <TNtuple.h>
  20 #include <TH1F.h>
  21 #include <TH2F.h>
  22 #include <TH3F.h>
  23 #include <TChain.h>
  24 #include "AliESDInputHandlerRP.h"
  25 #include "AliAnalysisTaskCheckHFMCProd.h"
  26
  27 /**************************************************************************
  28  * Copyright(c) 1998-2012, ALICE Experiment at CERN, All rights reserved. *
  29  *                                                                        *
  30  * Author: The ALICE Off-line Project.                                    *
  31  * Contributors are mentioned in the code where appropriate.              *
  32  *                                                                        *
  33  * Permission to use, copy, modify and distribute this software and its   *
  34  * documentation strictly for non-commercial purposes is hereby granted   *
  35  * without fee, provided that the above copyright notice appears in all   *
  36  * copies and that both the copyright notice and this permission notice   *
  37  * appear in the supporting documentation. The authors make no claims     *
  38  * about the suitability of this software for any purpose. It is          *
  39  * provided "as is" without express or implied warranty.                  *
  40  **************************************************************************/
  41
  42 /* $Id$ */
  43
  44 //*************************************************************************
  45 // Implementation of class AliAnalysisTaskCheckHFMCProd
  46 // AliAnalysisTask to check MC production at ESD+Kine level
  47 //
  48 //
  49 // Authors: F. Prino, prino@to.infn.it
  50 //
  51 //*************************************************************************
  52
  53 ClassImp(AliAnalysisTaskCheckHFMCProd)
  54 //______________________________________________________________________________
  55 AliAnalysisTaskCheckHFMCProd::AliAnalysisTaskCheckHFMCProd() : AliAnalysisTaskSE("HFMCChecks"),
  56   fOutput(0),
  57   fHistoNEvents(0),
  58   fHistoPhysPrim(0),
  59   fHistoTracks(0),
  60   fHistoSelTracks(0),
  61   fHistoTracklets(0),
  62   fHistoTrackletsEta1(0),
  63   fHistoPtPhysPrim(0),
  64   fHistoEtaPhysPrim(0),
  65   fHistoSPD3DVtxX(0),
  66   fHistoSPD3DVtxY(0),
  67   fHistoSPD3DVtxZ(0),
  68   fHistoSPDZVtxX(0),
  69   fHistoSPDZVtxY(0),
  70   fHistoSPDZVtxZ(0),
  71   fHistoTRKVtxX(0),
  72   fHistoTRKVtxY(0),
  73   fHistoTRKVtxZ(0),
  74   fHistoNcharmed(0),
  75   fHistoNbVsNc(0),
  76   fHistOriginPrompt(0),
  77   fHistOriginFeeddown(0),
  78   fHistMotherID(0),
  79   fHistDSpecies(0),
  80   fHistBSpecies(0),
  81   fHistNcollHFtype(0),
  82   fHistEtaPhiPtGenEle(0),
  83   fHistEtaPhiPtGenPi(0),
  84   fHistEtaPhiPtGenK(0),
  85   fHistEtaPhiPtGenPro(0),
  86   fHistEtaPhiPtRecEle(0),
  87   fHistEtaPhiPtRecPi(0),
  88   fHistEtaPhiPtRecK(0),
  89   fHistEtaPhiPtRecPro(0),
  90   fSearchUpToQuark(kFALSE),
  91   fSystem(0),
  92   fESDtrackCuts(0x0),
  93   fReadMC(kTRUE)
  94 {
  95   //
  96   DefineInput(0, TChain::Class());
  97   DefineOutput(1, TList::Class());
  98 }
  99
 100
 101 //___________________________________________________________________________
 102 AliAnalysisTaskCheckHFMCProd::~AliAnalysisTaskCheckHFMCProd(){
 103   //
 104   if (AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return;
 105   if (fOutput) {
 106     delete fOutput;
 107     fOutput = 0;
 108   }
 109   delete fESDtrackCuts;
 110
 111 }
 112
 113 //___________________________________________________________________________
 114 void AliAnalysisTaskCheckHFMCProd::UserCreateOutputObjects() {
 115   // create output histos
 116
 117   fOutput = new TList();
 118   fOutput->SetOwner();
 119   fOutput->SetName("OutputHistos");
 120
 121   fHistoNEvents = new TH1F("hNEvents", "Number of processed events",3,-0.5,2.5);
 122   fHistoNEvents->Sumw2();
 123   fHistoNEvents->SetMinimum(0);
 124   fOutput->Add(fHistoNEvents);
 125
 126   Double_t maxMult=100.;
 127   if(fSystem==1) maxMult=10000.;
 128   if(fSystem==2) maxMult=500.;
 129   fHistoPhysPrim = new TH1F("hPhysPrim","",100,-0.5,maxMult-0.5);
 130   fHistoPhysPrim->Sumw2();
 131   fOutput->Add(fHistoPhysPrim);
 132   fHistoTracks = new TH1F("hTracks","",100,-0.5,maxMult*2-0.5);
 133   fHistoTracks->Sumw2();
 134   fOutput->Add(fHistoTracks);
 135   fHistoSelTracks = new TH1F("hSelTracks","",100,-0.5,maxMult-0.5);
 136   fHistoSelTracks->Sumw2();
 137   fOutput->Add(fHistoSelTracks);
 138   fHistoTracklets = new TH1F("hTracklets","",100,-0.5,maxMult-0.5);
 139   fHistoTracklets->Sumw2();
 140   fOutput->Add(fHistoTracklets);
 141   fHistoTrackletsEta1 = new TH1F("hTrackletsEta1","",100,-0.5,maxMult-0.5);
 142   fHistoTrackletsEta1->Sumw2();
 143   fOutput->Add(fHistoTrackletsEta1);
 144   fHistoPtPhysPrim = new TH1F("hPtPhysPrim","",100,0.,20.);
 145   fHistoPtPhysPrim->Sumw2();
 146   fOutput->Add(fHistoPtPhysPrim);
 147   fHistoEtaPhysPrim = new TH1F("hEtaPhysPrim","",100,-10.,10.);
 148   fHistoEtaPhysPrim->Sumw2();
 149   fOutput->Add(fHistoEtaPhysPrim);
 150
 151   fHistoSPD3DVtxX = new TH1F("hSPD3DvX","",100,-1.,1.);
 152   fHistoSPD3DVtxX->Sumw2();
 153   fOutput->Add(fHistoSPD3DVtxX);
 154   fHistoSPD3DVtxY = new TH1F("hSPD3DvY","",100,-1.,1.);
 155   fHistoSPD3DVtxY->Sumw2();
 156   fOutput->Add(fHistoSPD3DVtxY);
 157   fHistoSPD3DVtxZ = new TH1F("hSPD3DvZ","",100,-15.,15.);
 158   fHistoSPD3DVtxZ->Sumw2();
 159   fOutput->Add(fHistoSPD3DVtxZ);
 160
 161   fHistoSPDZVtxX = new TH1F("hSPDZvX","",100,-1.,1.);
 162   fHistoSPDZVtxX->Sumw2();
 163   fOutput->Add(fHistoSPDZVtxX);
 164   fHistoSPDZVtxY = new TH1F("hSPDZvY","",100,-1.,1.);
 165   fHistoSPDZVtxY->Sumw2();
 166   fOutput->Add(fHistoSPDZVtxY);
 167   fHistoSPDZVtxZ = new TH1F("hSPDZvZ","",100,-15.,15.);
 168   fHistoSPDZVtxZ->Sumw2();
 169   fOutput->Add(fHistoSPDZVtxZ);
 170
 171
 172   fHistoTRKVtxX = new TH1F("hTRKvX","",100,-1.,1.);
 173   fHistoTRKVtxX->Sumw2();
 174   fOutput->Add(fHistoTRKVtxX);
 175   fHistoTRKVtxY = new TH1F("hTRKvY","",100,-1.,1.);
 176   fHistoTRKVtxY->Sumw2();
 177   fOutput->Add(fHistoTRKVtxY);
 178   fHistoTRKVtxZ = new TH1F("hTRKvZ","",100,-15.,15.);
 179   fHistoTRKVtxZ->Sumw2();
 180   fOutput->Add(fHistoTRKVtxZ);
 181
 182   Int_t nBinscb=11;
 183   if(fSystem==1) nBinscb=200;
 184   if(fSystem==2) nBinscb=21;
 185   Double_t maxncn=nBinscb-0.5;
 186   fHistoNcharmed = new TH2F("hncharmed","",100,-0.5,maxMult-0.5,nBinscb,-0.5,maxncn);
 187   fHistoNcharmed->Sumw2();
 188   fOutput->Add(fHistoNcharmed);
 189   fHistoNbVsNc = new TH2F("hnbvsnc","",nBinscb,-0.5,maxncn,nBinscb,-0.5,maxncn);
 190   fHistoNbVsNc->Sumw2();
 191   fOutput->Add(fHistoNbVsNc);
 192
 193   fHistYPtPrompt[0] = new TH2F("hyptD0prompt","D0 - Prompt",40,0.,40.,20,-2.,2.);
 194   fHistYPtPrompt[1] = new TH2F("hyptDplusprompt","Dplus - Prompt",40,0.,40.,20,-2.,2.);
 195   fHistYPtPrompt[2] = new TH2F("hyptDstarprompt","Dstar - Prompt",40,0.,40.,20,-2.,2.);
 196   fHistYPtPrompt[3] = new TH2F("hyptDsprompt","Ds - Prompt",40,0.,40.,20,-2.,2.);
 197   fHistYPtPrompt[4] = new TH2F("hyptLcprompt","Lc - Prompt",40,0.,40.,20,-2.,2.);
 198
 199   fHistBYPtAllDecay[0] = new TH2F("hyptB0AllDecay","B0 - All",40,0.,40.,40,-2.,2.);
 200   fHistBYPtAllDecay[1] = new TH2F("hyptBplusAllDecay","Bplus - All",40,0.,40.,40,-2.,2.);
 201   fHistBYPtAllDecay[2] = new TH2F("hyptBstarAllDecay","Bstar - All",40,0.,40.,40,-2.,2.);
 202   fHistBYPtAllDecay[3] = new TH2F("hyptBsAllDecay","Bs - All",40,0.,40.,40,-2.,2.);
 203   fHistBYPtAllDecay[4] = new TH2F("hyptLbAllDecay","LB - All",40,0.,40.,40,-2.,2.);
 204
 205   fHistYPtAllDecay[0] = new TH2F("hyptD0AllDecay","D0 - All",40,0.,40.,40,-2.,2.);
 206   fHistYPtAllDecay[1] = new TH2F("hyptDplusAllDecay","Dplus - All",40,0.,40.,40,-2.,2.);
 207   fHistYPtAllDecay[2] = new TH2F("hyptDstarAllDecay","Dstar - All",40,0.,40.,40,-2.,2.);
 208   fHistYPtAllDecay[3] = new TH2F("hyptDsAllDecay","Ds - All",40,0.,40.,40,-2.,2.);
 209   fHistYPtAllDecay[4] = new TH2F("hyptLcAllDecay","Lc - All",40,0.,40.,40,-2.,2.);
 210
 211   fHistYPtPromptAllDecay[0] = new TH2F("hyptD0promptAllDecay","D0 - Prompt",40,0.,40.,40,-2.,2.);
 212   fHistYPtPromptAllDecay[1] = new TH2F("hyptDpluspromptAllDecay","Dplus - Prompt",40,0.,40.,40,-2.,2.);
 213   fHistYPtPromptAllDecay[2] = new TH2F("hyptDstarpromptAllDecay","Dstar - Prompt",40,0.,40.,40,-2.,2.);
 214   fHistYPtPromptAllDecay[3] = new TH2F("hyptDspromptAllDecay","Ds - Prompt",40,0.,40.,40,-2.,2.);
 215   fHistYPtPromptAllDecay[4] = new TH2F("hyptLcpromptAllDecay","Lc - Prompt",40,0.,40.,40,-2.,2.);
 216
 217   fHistYPtFeeddownAllDecay[0] = new TH2F("hyptD0feeddownAllDecay","D0 - FromB",40,0.,40.,40,-2.,2.);
 218   fHistYPtFeeddownAllDecay[1] = new TH2F("hyptDplusfeeddownAllDecay","Dplus - FromB",40,0.,40.,40,-2.,2.);
 219   fHistYPtFeeddownAllDecay[2] = new TH2F("hyptDstarfeeddownAllDecay","Dstar - FromB",40,0.,40.,40,-2.,2.);
 220   fHistYPtFeeddownAllDecay[3] = new TH2F("hyptDsfeeddownAllDecay","Ds - FromB",40,0.,40.,40,-2.,2.);
 221   fHistYPtFeeddownAllDecay[4] = new TH2F("hyptLcfeeddownAllDecay","Lc - FromB",40,0.,40.,40,-2.,2.);
 222
 223
 224   fHistYPtFeeddown[0] = new TH2F("hyptD0feeddown","D0 - Feeddown",40,0.,40.,20,-2.,2.);
 225   fHistYPtFeeddown[1] = new TH2F("hyptDplusfeeddown","Dplus - Feeddown",40,0.,40.,20,-2.,2.);
 226   fHistYPtFeeddown[2] = new TH2F("hyptDstarfeedown","Dstar - Feeddown",40,0.,40.,20,-2.,2.);
 227   fHistYPtFeeddown[3] = new TH2F("hyptDsfeedown","Ds - Feeddown",40,0.,40.,20,-2.,2.);
 228   fHistYPtFeeddown[4] = new TH2F("hyptLcfeedown","Lc - Feeddown",40,0.,40.,20,-2.,2.);
 229
 230   for(Int_t ih=0; ih<5; ih++){
 231     fHistBYPtAllDecay[ih]->Sumw2();
 232     fHistBYPtAllDecay[ih]->SetMinimum(0);
 233     fOutput->Add(fHistBYPtAllDecay[ih]);
 234     fHistYPtAllDecay[ih]->Sumw2();
 235     fHistYPtAllDecay[ih]->SetMinimum(0);
 236     fOutput->Add(fHistYPtAllDecay[ih]);
 237     fHistYPtPromptAllDecay[ih]->Sumw2();
 238     fHistYPtPromptAllDecay[ih]->SetMinimum(0);
 239     fOutput->Add(fHistYPtPromptAllDecay[ih]);
 240     fHistYPtFeeddownAllDecay[ih]->Sumw2();
 241     fHistYPtFeeddownAllDecay[ih]->SetMinimum(0);
 242     fOutput->Add(fHistYPtFeeddownAllDecay[ih]);
 243     fHistYPtPrompt[ih]->Sumw2();
 244     fHistYPtPrompt[ih]->SetMinimum(0);
 245     fOutput->Add(fHistYPtPrompt[ih]);
 246     fHistYPtFeeddown[ih]->Sumw2();
 247     fHistYPtFeeddown[ih]->SetMinimum(0);
 248     fOutput->Add(fHistYPtFeeddown[ih]);
 249   }
 250
 251   fHistYPtD0byDecChannel[0] = new TH2F("hyptD02","D0 - 2prong",40,0.,40.,20,-2.,2.);
 252   fHistYPtD0byDecChannel[1] = new TH2F("hyptD04","D0 - 4prong",40,0.,40.,20,-2.,2.);
 253   fHistYPtDplusbyDecChannel[0] = new TH2F("hyptDplusnonreson","Dplus - non reson",40,0.,40.,20,-2.,2.);
 254   fHistYPtDplusbyDecChannel[1] = new TH2F("hyptDplusreson","Dplus - reson via K0*",40,0.,40.,20,-2.,2.);
 255   fHistYPtDsbyDecChannel[0] = new TH2F("hyptDsphi","Ds - vis Phi",40,0.,40.,20,-2.,2.);
 256   fHistYPtDsbyDecChannel[1] = new TH2F("hyptDsk0st","Ds - via k0*",40,0.,40.,20,-2.,2.);
 257
 258   for(Int_t ih=0; ih<2; ih++){
 259
 260     fHistYPtD0byDecChannel[ih]->Sumw2();
 261     fHistYPtD0byDecChannel[ih]->SetMinimum(0);
 262     fOutput->Add(fHistYPtD0byDecChannel[ih]);
 263     fHistYPtDplusbyDecChannel[ih]->Sumw2();
 264     fHistYPtDplusbyDecChannel[ih]->SetMinimum(0);
 265     fOutput->Add(fHistYPtDplusbyDecChannel[ih]);
 266     fHistYPtDsbyDecChannel[ih]->Sumw2();
 267     fHistYPtDsbyDecChannel[ih]->SetMinimum(0);
 268     fOutput->Add(fHistYPtDsbyDecChannel[ih]);
 269   }
 270
 271   fHistOriginPrompt=new TH1F("hOriginPrompt","",100,0.,0.5);
 272   fHistOriginPrompt->Sumw2();
 273   fHistOriginPrompt->SetMinimum(0);
 274   fOutput->Add(fHistOriginPrompt);
 275   fHistOriginFeeddown=new TH1F("hOriginFeeddown","",100,0.,0.5);
 276   fHistOriginFeeddown->Sumw2();
 277   fHistOriginFeeddown->SetMinimum(0);
 278   fOutput->Add(fHistOriginFeeddown);
 279   fHistMotherID=new TH1F("hMotherID","",1000,-1.5,998.5);
 280   fHistMotherID->SetMinimum(0);
 281   fOutput->Add(fHistMotherID);
 282   fHistDSpecies=new TH1F("hDSpecies","",10,-0.5,9.5);
 283   fHistDSpecies->GetXaxis()->SetBinLabel(1,"D0");
 284   fHistDSpecies->GetXaxis()->SetBinLabel(2,"D0bar");
 285   fHistDSpecies->GetXaxis()->SetBinLabel(3,"D+");
 286   fHistDSpecies->GetXaxis()->SetBinLabel(4,"D-");
 287   fHistDSpecies->GetXaxis()->SetBinLabel(5,"D*+");
 288   fHistDSpecies->GetXaxis()->SetBinLabel(6,"D*-");
 289   fHistDSpecies->GetXaxis()->SetBinLabel(7,"Ds+");
 290   fHistDSpecies->GetXaxis()->SetBinLabel(8,"Ds-");
 291   fHistDSpecies->GetXaxis()->SetBinLabel(9,"Lc+");
 292   fHistDSpecies->GetXaxis()->SetBinLabel(10,"Lc-");
 293   fHistDSpecies->SetMinimum(0);
 294   fOutput->Add(fHistDSpecies);
 295   fHistBSpecies=new TH1F("hBSpecies","",10,-0.5,9.5);
 296   fHistBSpecies->GetXaxis()->SetBinLabel(1,"B0");
 297   fHistBSpecies->GetXaxis()->SetBinLabel(2,"B0bar");
 298   fHistBSpecies->GetXaxis()->SetBinLabel(3,"B+");
 299   fHistBSpecies->GetXaxis()->SetBinLabel(4,"B-");
 300   fHistBSpecies->GetXaxis()->SetBinLabel(5,"B*+");
 301   fHistBSpecies->GetXaxis()->SetBinLabel(6,"B*-");
 302   fHistBSpecies->GetXaxis()->SetBinLabel(7,"Bs+");
 303   fHistBSpecies->GetXaxis()->SetBinLabel(8,"Bs-");
 304   fHistBSpecies->GetXaxis()->SetBinLabel(9,"Lb+");
 305   fHistBSpecies->GetXaxis()->SetBinLabel(10,"Lb-");
 306   fHistBSpecies->SetMinimum(0);
 307   fOutput->Add(fHistBSpecies);
 308
 309   fHistNcollHFtype=new TH2F("hNcollHFtype","",5,-1.5,3.5,30,-0.5,29.5);
 310   fOutput->Add(fHistNcollHFtype);
 311
 312   Double_t binseta[11]={-1.0,-0.8,-0.6,-0.4,-0.2,0.,0.2,0.4,0.6,0.8,1.0};
 313   const Int_t nBinsPhi=40;
 314   Double_t binsphi[nBinsPhi+1];
 315   for(Int_t ib=0; ib<=nBinsPhi; ib++) binsphi[ib]=ib*TMath::Pi()/20.;
 316   const Int_t nBinsPt=24;
 317   Double_t binspt[nBinsPt+1]={0.,0.10,0.15,0.2,0.25,
 318                               0.3,0.4,0.5,0.6,0.7,
 319                               0.8,0.9,1.,1.25,1.5,
 320                               1.75,2.,2.5,3.,4.,
 321                               5.,7.5,10.,15.,20.};
 322
 323    fHistEtaPhiPtGenEle=new TH3F("hEtaPhiPtGenEle","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
 324   fOutput->Add(fHistEtaPhiPtGenEle);
 325   fHistEtaPhiPtGenPi=new TH3F("hEtaPhiPtGenPi","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
 326   fOutput->Add(fHistEtaPhiPtGenPi);
 327   fHistEtaPhiPtGenK=new TH3F("hEtaPhiPtGenK","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
 328   fOutput->Add(fHistEtaPhiPtGenK);
 329   fHistEtaPhiPtGenPro=new TH3F("hEtaPhiPtGenPro","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
 330   fOutput->Add(fHistEtaPhiPtGenPro);
 331
 332
 333   fHistEtaPhiPtRecEle=new TH3F("hEtaPhiPtRecEle","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
 334   fOutput->Add(fHistEtaPhiPtRecEle);
 335   fHistEtaPhiPtRecPi=new TH3F("hEtaPhiPtRecPi","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
 336   fOutput->Add(fHistEtaPhiPtRecPi);
 337   fHistEtaPhiPtRecK=new TH3F("hEtaPhiPtRecK","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
 338   fOutput->Add(fHistEtaPhiPtRecK);
 339   fHistEtaPhiPtRecPro=new TH3F("hEtaPhiPtRecPro","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
 340   fOutput->Add(fHistEtaPhiPtRecPro);
 341
 342
 343
 344   PostData(1,fOutput);
 345
 346 }
 347 //______________________________________________________________________________
 348 void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *)
 349 {
 350   //
 351
 352   AliESDEvent *esd = (AliESDEvent*) (InputEvent());
 353
 354
 355   if(!esd) {
 356     printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
 357     return;
 358   }
 359
 360   fHistoNEvents->Fill(0);
 361
 362   if(!fESDtrackCuts){
 363     Int_t year=2011;
 364     if(esd->GetRunNumber()<=139517) year=2010;
 365     if(year==2010) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kFALSE);
 366     else fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
 367     fESDtrackCuts->SetMaxDCAToVertexXY(2.4);
 368     fESDtrackCuts->SetMaxDCAToVertexZ(3.2);
 369     fESDtrackCuts->SetDCAToVertex2D(kTRUE);
 370     fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 371                                             AliESDtrackCuts::kAny);
 372   }
 373
 374   Int_t nTracks=esd->GetNumberOfTracks();
 375   fHistoTracks->Fill(nTracks);
 376   Int_t nSelTracks=0;
 377   for(Int_t it=0; it<nTracks; it++){
 378     AliESDtrack* tr=esd->GetTrack(it);
 379     UInt_t status=tr->GetStatus();
 380     if(!(status&AliESDtrack::kITSrefit)) continue;
 381     if(!(status&AliESDtrack::kTPCin)) continue;
 382     nSelTracks++;
 383   }
 384   fHistoSelTracks->Fill(nSelTracks);
 385
 386   const AliMultiplicity* mult=esd->GetMultiplicity();
 387   Int_t nTracklets=mult->GetNumberOfTracklets();
 388   Int_t nTracklets1=0;
 389   for(Int_t it=0; it<nTracklets; it++){
 390     Double_t eta=TMath::Abs(mult->GetEta(it));
 391     if(eta<1) nTracklets1++;
 392   }
 393   fHistoTracklets->Fill(nTracklets);
 394   fHistoTrackletsEta1->Fill(nTracklets1);
 395
 396   const AliESDVertex *spdv=esd->GetVertex();
 397   if(spdv && spdv->IsFromVertexer3D()){
 398     fHistoSPD3DVtxX->Fill(spdv->GetX());
 399     fHistoSPD3DVtxY->Fill(spdv->GetY());
 400     fHistoSPD3DVtxZ->Fill(spdv->GetZ());
 401   }
 402   if(spdv && spdv->IsFromVertexerZ()){
 403     fHistoSPDZVtxX->Fill(spdv->GetX());
 404     fHistoSPDZVtxY->Fill(spdv->GetY());
 405     fHistoSPDZVtxZ->Fill(spdv->GetZ());
 406   }
 407   const AliESDVertex *trkv=esd->GetPrimaryVertex();
 408   if(trkv && trkv->GetNContributors()>1){
 409     fHistoTRKVtxX->Fill(trkv->GetX());
 410     fHistoTRKVtxY->Fill(trkv->GetY());
 411     fHistoTRKVtxZ->Fill(trkv->GetZ());
 412   }
 413
 414   AliStack* stack=0;
 415   if(fReadMC){
 416     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 417     if (!eventHandler) {
 418       Printf("ERROR: Could not retrieve MC event handler");
 419       return;
 420     }
 421     AliMCEvent* mcEvent = eventHandler->MCEvent();
 422     if (!mcEvent) {
 423       Printf("ERROR: Could not retrieve MC event");
 424       return;
 425     }
 426     stack = mcEvent->Stack();
 427     if (!stack) {
 428       Printf("ERROR: stack not available");
 429       return;
 430     }
 431     const AliVVertex* mcVert=mcEvent->GetPrimaryVertex();
 432     if(!mcVert){
 433       Printf("ERROR: generated vertex not available");
 434       return;
 435     }
 436     if(TMath::Abs(mcVert->GetZ())>10) return;
 437
 438     //    const AliHeader* h=(AliHeader*)mcEvent->GetHeader();
 439     //    cout<<h<<endl;
 440     TString genname=mcEvent->GenEventHeader()->ClassName();
 441     Int_t nColl=-1;
 442     Int_t typeHF=-1;
 443     TList* lgen=0x0;
 444     if(genname.Contains("CocktailEventHeader")){
 445       AliGenCocktailEventHeader *cockhead=(AliGenCocktailEventHeader*)mcEvent->GenEventHeader();
 446       lgen=cockhead->GetHeaders();
 447       for(Int_t ig=0; ig<lgen->GetEntries(); ig++){
 448         AliGenerator* gen=(AliGenerator*)lgen->At(ig);
 449         TString title=gen->GetName();
 450         if(title.Contains("bchadr")) typeHF=1;
 451         else if(title.Contains("chadr")) typeHF=0;
 452         else if(title.Contains("bele")) typeHF=3;
 453         else if(title.Contains("cele")) typeHF=2;
 454       }
 455       nColl=lgen->GetEntries();
 456       fHistNcollHFtype->Fill(typeHF,nColl);
 457     }
 458     Int_t nParticles=stack->GetNtrack();
 459     Double_t dNchdy = 0.;
 460     Int_t nb = 0, nc=0;
 461     Int_t nCharmed=0;
 462     Int_t nPhysPrim=0;
 463     for (Int_t i=0;i<nParticles;i++){
 464       TParticle* part = (TParticle*)stack->Particle(i);
 465       Int_t absPdg=TMath::Abs(part->GetPdgCode());
 466       Int_t pdg=part->GetPdgCode();
 467       if(absPdg==4) nc++;
 468       if(absPdg==5) nb++;
 469       if(stack->IsPhysicalPrimary(i)){
 470         Double_t eta=part->Eta();
 471         fHistoEtaPhysPrim->Fill(eta);
 472         if(absPdg==11) fHistEtaPhiPtGenEle->Fill(eta,part->Phi(),part->Pt());
 473         else if(absPdg==211) fHistEtaPhiPtGenPi->Fill(eta,part->Phi(),part->Pt());
 474         else if(absPdg==321) fHistEtaPhiPtGenK->Fill(eta,part->Phi(),part->Pt());
 475         else if(absPdg==2212) fHistEtaPhiPtGenPro->Fill(eta,part->Phi(),part->Pt());
 476
 477         if(TMath::Abs(eta)<0.5){
 478           dNchdy+=0.6666;   // 2/3 for the ratio charged/all
 479           nPhysPrim++;
 480         }
 481         if(TMath::Abs(eta)<0.9){
 482           fHistoPtPhysPrim->Fill(part->Pt());
 483         }
 484       }
 485       Float_t rapid=-999.;
 486       if (part->Energy() != TMath::Abs(part->Pz())){
 487         rapid=0.5*TMath::Log((part->Energy()+part->Pz())/(part->Energy()-part->Pz()));
 488       }
 489
 490
 491       Int_t iPart=-1;
 492       Int_t iType=0;
 493       Int_t iSpecies=-1;
 494       Int_t dummy[4];
 495       if(absPdg==421){
 496         iSpecies=0;
 497         iType=AliVertexingHFUtils::CheckD0Decay(stack,i,dummy);
 498         if(iType>0) iPart=0;
 499       }
 500       else if(absPdg==411){
 501         iSpecies=1;
 502         iType=AliVertexingHFUtils::CheckDplusDecay(stack,i,dummy);
 503         if(iType>0) iPart=1;
 504       }
 505       else if(absPdg==413){
 506         iSpecies=2;
 507         iType=AliVertexingHFUtils::CheckDstarDecay(stack,i,dummy);
 508         if(iType>0) iPart=2;
 509       }
 510       else if(absPdg==431){
 511         iSpecies=3;
 512         iType=AliVertexingHFUtils::CheckDsDecay(stack,i,dummy);
 513         if(iType==1 || iType==2) iPart=3;
 514       }
 515       else if(absPdg==4122){
 516         iSpecies=4;
 517         iType=CheckLcDecay(i,stack);
 518         if(iType>=0) iPart=4;
 519       }
 520       if(iSpecies>=0) fHistYPtAllDecay[iSpecies]->Fill(part->Pt(),rapid);
 521
 522       // check beauty mesons
 523       if(absPdg==511) fHistBYPtAllDecay[0]->Fill(part->Pt(),rapid);
 524       else if(absPdg==521) fHistBYPtAllDecay[1]->Fill(part->Pt(),rapid);
 525       else if(absPdg==513) fHistBYPtAllDecay[2]->Fill(part->Pt(),rapid);
 526       else if(absPdg==531) fHistBYPtAllDecay[3]->Fill(part->Pt(),rapid);
 527       else if(absPdg==5122) fHistBYPtAllDecay[4]->Fill(part->Pt(),rapid);
 528
 529       if(pdg==511) fHistBSpecies->Fill(0);
 530       else if(pdg==-511) fHistBSpecies->Fill(1);
 531       else if(pdg==521) fHistBSpecies->Fill(2);
 532       else if(pdg==-521) fHistBSpecies->Fill(3);
 533       else if(pdg==513) fHistBSpecies->Fill(4);
 534       else if(pdg==-513) fHistBSpecies->Fill(5);
 535       else if(pdg==531) fHistBSpecies->Fill(6);
 536       else if(pdg==-531) fHistBSpecies->Fill(7);
 537       else if(pdg==5122) fHistBSpecies->Fill(8);
 538       else if(pdg==-5122) fHistBSpecies->Fill(9);
 539
 540      if(iSpecies<0) continue; // not a charm meson
 541
 542       if(pdg==421) fHistDSpecies->Fill(0);
 543       else if(pdg==-421) fHistDSpecies->Fill(1);
 544       else if(pdg==411) fHistDSpecies->Fill(2);
 545       else if(pdg==-411) fHistDSpecies->Fill(3);
 546       else if(pdg==413) fHistDSpecies->Fill(4);
 547       else if(pdg==-413) fHistDSpecies->Fill(5);
 548       else if(pdg==431) fHistDSpecies->Fill(6);
 549       else if(pdg==-431) fHistDSpecies->Fill(7);
 550       else if(pdg==4122) fHistDSpecies->Fill(8);
 551       else if(pdg==-4122) fHistDSpecies->Fill(9);
 552
 553       Double_t distx=part->Vx()-mcVert->GetX();
 554       Double_t disty=part->Vy()-mcVert->GetY();
 555       Double_t distz=part->Vz()-mcVert->GetZ();
 556       Double_t distToVert=TMath::Sqrt(distx*distx+disty*disty+distz*distz);
 557       fHistMotherID->Fill(part->GetFirstMother());
 558       Int_t iFromB=AliVertexingHFUtils::CheckOrigin(stack,part,fSearchUpToQuark);
 559       if(iFromB==4){
 560         fHistYPtPromptAllDecay[iSpecies]->Fill(part->Pt(),rapid);
 561         fHistOriginPrompt->Fill(distToVert);
 562       }
 563       else if(iFromB==5){
 564         fHistYPtFeeddownAllDecay[iSpecies]->Fill(part->Pt(),rapid);
 565         fHistOriginFeeddown->Fill(distToVert);
 566       }
 567
 568       if(iPart<0) continue;
 569       if(iType<0) continue;
 570       nCharmed++;
 571       if(iPart==0 && iType>0 && iType<=2){
 572         fHistYPtD0byDecChannel[iType-1]->Fill(part->Pt(),rapid);
 573       }else if(iPart==1 && iType>0 && iType<=2){
 574         fHistYPtDplusbyDecChannel[iType-1]->Fill(part->Pt(),rapid);
 575       }else if(iPart==3 &&  iType>0 && iType<=2){
 576         fHistYPtDsbyDecChannel[iType-1]->Fill(part->Pt(),rapid);
 577       }
 578
 579       if(iFromB==4 && iPart>=0 && iPart<5) fHistYPtPrompt[iPart]->Fill(part->Pt(),rapid);
 580       else if(iFromB==5 && iPart>=0 && iPart<5) fHistYPtFeeddown[iPart]->Fill(part->Pt(),rapid);
 581     }
 582
 583     for(Int_t i=0; i<nTracks; i++){
 584       AliESDtrack* track=esd->GetTrack(i);
 585       if(fESDtrackCuts->AcceptTrack(track)){
 586         Int_t label=TMath::Abs(track->GetLabel());
 587         if(stack->IsPhysicalPrimary(label)){
 588           TParticle* part = (TParticle*)stack->Particle(label);
 589           Int_t absPdg=TMath::Abs(part->GetPdgCode());
 590           Double_t eta=part->Eta();
 591           if(absPdg==11) fHistEtaPhiPtRecEle->Fill(eta,part->Phi(),part->Pt());
 592           else if(absPdg==211) fHistEtaPhiPtRecPi->Fill(eta,part->Phi(),part->Pt());
 593           else if(absPdg==321) fHistEtaPhiPtRecK->Fill(eta,part->Phi(),part->Pt());
 594           else if(absPdg==2212) fHistEtaPhiPtRecPro->Fill(eta,part->Phi(),part->Pt());
 595         }
 596       }
 597     }
 598     fHistoNcharmed->Fill(dNchdy,nCharmed);
 599     fHistoNbVsNc->Fill(nc,nb);
 600     fHistoPhysPrim->Fill(nPhysPrim);
 601   }
 602
 603   PostData(1,fOutput);
 604
 605 }
 606 //______________________________________________________________________________
 607 void AliAnalysisTaskCheckHFMCProd::Terminate(Option_t */*option*/)
 608 {
 609   // Terminate analysis
 610   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 611   if (!fOutput) {
 612     printf("ERROR: fOutput not available\n");
 613     return;
 614   }
 615
 616   return;
 617 }
 618
 619
 620
 621
 622 //______________________________________________________________________________
 623 Int_t AliAnalysisTaskCheckHFMCProd::CheckLcDecay(Int_t labLc, AliStack* stack) const{
 624   if(labLc<0) return -1;
 625   TParticle* dp = (TParticle*)stack->Particle(labLc);
 626   Int_t pdgdp=dp->GetPdgCode();
 627   Int_t nDau=dp->GetNDaughters();
 628
 629   if(nDau==3){
 630     Int_t nKaons=0;
 631     Int_t nPions=0;
 632     Int_t nProtons=0;
 633     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 634       if(iDau<0) return -1;
 635       TParticle* dau=(TParticle*)stack->Particle(iDau);
 636       Int_t pdgdau=dau->GetPdgCode();
 637       if(TMath::Abs(pdgdau)==321){
 638         if(pdgdp>0 && pdgdau>0) return -1;
 639         if(pdgdp<0 && pdgdau<0) return -1;
 640         nKaons++;
 641       }else if(TMath::Abs(pdgdau)==211){
 642         if(pdgdp<0 && pdgdau>0) return -1;
 643         if(pdgdp>0 && pdgdau<0) return -1;
 644         nPions++;
 645       }else if(TMath::Abs(pdgdau)==2212){
 646         if(pdgdp<0 && pdgdau>0) return -1;
 647         if(pdgdp>0 && pdgdau<0) return -1;
 648         nProtons++;
 649       }else{
 650         return -1;
 651       }
 652     }
 653     if(nPions!=1) return -1;
 654     if(nKaons!=1) return -1;
 655     if(nProtons!=1) return -1;
 656     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 657       if(iDau<0) return -1;
 658     }
 659     return 0;
 660   }
 661
 662   if(nDau==2){
 663     Int_t nKaons=0;
 664     Int_t nPions=0;
 665     Int_t nProtons=0;
 666     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 667       if(iDau<0) return -1;
 668       TParticle* dau=(TParticle*)stack->Particle(iDau);
 669       Int_t pdgdau=dau->GetPdgCode();
 670       if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 || TMath::Abs(pdgdau)==2224
 671          || TMath::Abs(pdgdau)==3122 || TMath::Abs(pdgdau)==311){
 672         Int_t nDauRes=dau->GetNDaughters();
 673         if(nDauRes!=2)  return -1;
 674         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 675           if(resDau<0) return -1;
 676           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 677           Int_t pdgresdau=resdau->GetPdgCode();
 678           if(TMath::Abs(pdgresdau)==321){
 679             if(pdgdp>0 && pdgresdau>0) return -1;
 680             if(pdgdp<0 && pdgresdau<0) return -1;
 681             nKaons++;
 682           }
 683           if(TMath::Abs(pdgresdau)==211){
 684             if(pdgdp<0 && pdgresdau>0) return -1;
 685             if(pdgdp>0 && pdgresdau<0) return -1;
 686             nPions++;
 687           }
 688           if(TMath::Abs(pdgresdau)==2212){
 689             if(pdgdp<0 && pdgresdau>0) return -1;
 690             if(pdgdp>0 && pdgresdau<0) return -1;
 691             nProtons++;
 692           }
 693         }
 694       }else if(TMath::Abs(pdgdau)==321){
 695         if(pdgdp>0 && pdgdau>0) return -1;
 696         if(pdgdp<0 && pdgdau<0) return -1;
 697         nKaons++;
 698       }else if(TMath::Abs(pdgdau)==211){
 699         if(pdgdp<0 && pdgdau>0) return -1;
 700         if(pdgdp>0 && pdgdau<0) return -1;
 701         nPions++;
 702       }else if(TMath::Abs(pdgdau)==2212){
 703         if(pdgdp<0 && pdgdau>0) return -1;
 704         if(pdgdp>0 && pdgdau<0) return -1;
 705         nProtons++;
 706       }else{
 707         return -1;
 708       }
 709     }
 710     if(nPions!=1) return -1;
 711     if(nKaons!=1) return -1;
 712     if(nProtons!=1) return -1;
 713     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 714       if(iDau<0) return -1;
 715     }
 716     return 1;
 717   }
 718   return -1;
 719 }
 720