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