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 "AliMultiplicity.h"
  15 #include <TParticle.h>
  16 #include <TSystem.h>
  17 #include <TTree.h>
  18 #include <TNtuple.h>
  19 #include <TH1F.h>
  20 #include <TH2F.h>
  21 #include <TH3F.h>
  22 #include <TChain.h>
  23 #include "AliESDInputHandlerRP.h"
  24 #include "AliAnalysisTaskCheckHFMCProd.h"
  25
  26 /**************************************************************************
  27  * Copyright(c) 1998-2012, ALICE Experiment at CERN, All rights reserved. *
  28  *                                                                        *
  29  * Author: The ALICE Off-line Project.                                    *
  30  * Contributors are mentioned in the code where appropriate.              *
  31  *                                                                        *
  32  * Permission to use, copy, modify and distribute this software and its   *
  33  * documentation strictly for non-commercial purposes is hereby granted   *
  34  * without fee, provided that the above copyright notice appears in all   *
  35  * copies and that both the copyright notice and this permission notice   *
  36  * appear in the supporting documentation. The authors make no claims     *
  37  * about the suitability of this software for any purpose. It is          *
  38  * provided "as is" without express or implied warranty.                  *
  39  **************************************************************************/
  40
  41 /* $Id$ */
  42
  43 //*************************************************************************
  44 // Implementation of class AliAnalysisTaskCheckHFMCProd
  45 // AliAnalysisTask to check MC production at ESD+Kine level
  46 //
  47 //
  48 // Authors: F. Prino, prino@to.infn.it
  49 //
  50 //*************************************************************************
  51
  52 ClassImp(AliAnalysisTaskCheckHFMCProd)
  53 //______________________________________________________________________________
  54 AliAnalysisTaskCheckHFMCProd::AliAnalysisTaskCheckHFMCProd() : AliAnalysisTaskSE("HFMCChecks"),
  55   fOutput(0),
  56   fHistoNEvents(0),
  57   fHistoPhysPrim(0),
  58   fHistoTracks(0),
  59   fHistoSelTracks(0),
  60   fHistoTracklets(0),
  61   fHistoTrackletsEta1(0),
  62   fHistoPtPhysPrim(0),
  63   fHistoEtaPhysPrim(0),
  64   fHistoSPD3DVtxX(0),
  65   fHistoSPD3DVtxY(0),
  66   fHistoSPD3DVtxZ(0),
  67   fHistoSPDZVtxX(0),
  68   fHistoSPDZVtxY(0),
  69   fHistoSPDZVtxZ(0),
  70   fHistoTRKVtxX(0),
  71   fHistoTRKVtxY(0),
  72   fHistoTRKVtxZ(0),
  73   fHistoNcharmed(0),
  74   fHistoNbVsNc(0),
  75   fHistOriginPrompt(0),
  76   fHistOriginFeeddown(0),
  77   fHistMotherID(0),
  78   fHistDSpecies(0),
  79   fHistBSpecies(0),
  80   fHistNcollHFtype(0),
  81   fHistEtaPhiPtGenEle(0),
  82   fHistEtaPhiPtGenPi(0),
  83   fHistEtaPhiPtGenK(0),
  84   fHistEtaPhiPtGenPro(0),
  85   fHistEtaPhiPtRecEle(0),
  86   fHistEtaPhiPtRecPi(0),
  87   fHistEtaPhiPtRecK(0),
  88   fHistEtaPhiPtRecPro(0),
  89   fSearchUpToQuark(kFALSE),
  90   fSystem(0),
  91   fESDtrackCuts(0x0),
  92   fReadMC(kTRUE)
  93 {
  94   //
  95   DefineInput(0, TChain::Class());
  96   DefineOutput(1, TList::Class());
  97 }
  98
  99
 100 //___________________________________________________________________________
 101 AliAnalysisTaskCheckHFMCProd::~AliAnalysisTaskCheckHFMCProd(){
 102   //
 103   if (AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return;
 104   if (fOutput) {
 105     delete fOutput;
 106     fOutput = 0;
 107   }
 108   delete fESDtrackCuts;
 109
 110 }
 111
 112 //___________________________________________________________________________
 113 void AliAnalysisTaskCheckHFMCProd::UserCreateOutputObjects() {
 114   // create output histos
 115
 116   fOutput = new TList();
 117   fOutput->SetOwner();
 118   fOutput->SetName("OutputHistos");
 119
 120   fHistoNEvents = new TH1F("hNEvents", "Number of processed events",3,-0.5,2.5);
 121   fHistoNEvents->Sumw2();
 122   fHistoNEvents->SetMinimum(0);
 123   fOutput->Add(fHistoNEvents);
 124
 125   Double_t maxMult=100.;
 126   if(fSystem==1) maxMult=10000.;
 127   if(fSystem==2) maxMult=500.;
 128   fHistoPhysPrim = new TH1F("hPhysPrim","",100,-0.5,maxMult-0.5);
 129   fHistoPhysPrim->Sumw2();
 130   fOutput->Add(fHistoPhysPrim);
 131   fHistoTracks = new TH1F("hTracks","",100,-0.5,maxMult*2-0.5);
 132   fHistoTracks->Sumw2();
 133   fOutput->Add(fHistoTracks);
 134   fHistoSelTracks = new TH1F("hSelTracks","",100,-0.5,maxMult-0.5);
 135   fHistoSelTracks->Sumw2();
 136   fOutput->Add(fHistoSelTracks);
 137   fHistoTracklets = new TH1F("hTracklets","",100,-0.5,maxMult-0.5);
 138   fHistoTracklets->Sumw2();
 139   fOutput->Add(fHistoTracklets);
 140   fHistoTrackletsEta1 = new TH1F("hTrackletsEta1","",100,-0.5,maxMult-0.5);
 141   fHistoTrackletsEta1->Sumw2();
 142   fOutput->Add(fHistoTrackletsEta1);
 143   fHistoPtPhysPrim = new TH1F("hPtPhysPrim","",100,0.,20.);
 144   fHistoPtPhysPrim->Sumw2();
 145   fOutput->Add(fHistoPtPhysPrim);
 146   fHistoEtaPhysPrim = new TH1F("hEtaPhysPrim","",100,-10.,10.);
 147   fHistoEtaPhysPrim->Sumw2();
 148   fOutput->Add(fHistoEtaPhysPrim);
 149
 150   fHistoSPD3DVtxX = new TH1F("hSPD3DvX","",100,-1.,1.);
 151   fHistoSPD3DVtxX->Sumw2();
 152   fOutput->Add(fHistoSPD3DVtxX);
 153   fHistoSPD3DVtxY = new TH1F("hSPD3DvY","",100,-1.,1.);
 154   fHistoSPD3DVtxY->Sumw2();
 155   fOutput->Add(fHistoSPD3DVtxY);
 156   fHistoSPD3DVtxZ = new TH1F("hSPD3DvZ","",100,-15.,15.);
 157   fHistoSPD3DVtxZ->Sumw2();
 158   fOutput->Add(fHistoSPD3DVtxZ);
 159
 160   fHistoSPDZVtxX = new TH1F("hSPDZvX","",100,-1.,1.);
 161   fHistoSPDZVtxX->Sumw2();
 162   fOutput->Add(fHistoSPDZVtxX);
 163   fHistoSPDZVtxY = new TH1F("hSPDZvY","",100,-1.,1.);
 164   fHistoSPDZVtxY->Sumw2();
 165   fOutput->Add(fHistoSPDZVtxY);
 166   fHistoSPDZVtxZ = new TH1F("hSPDZvZ","",100,-15.,15.);
 167   fHistoSPDZVtxZ->Sumw2();
 168   fOutput->Add(fHistoSPDZVtxZ);
 169
 170
 171   fHistoTRKVtxX = new TH1F("hTRKvX","",100,-1.,1.);
 172   fHistoTRKVtxX->Sumw2();
 173   fOutput->Add(fHistoTRKVtxX);
 174   fHistoTRKVtxY = new TH1F("hTRKvY","",100,-1.,1.);
 175   fHistoTRKVtxY->Sumw2();
 176   fOutput->Add(fHistoTRKVtxY);
 177   fHistoTRKVtxZ = new TH1F("hTRKvZ","",100,-15.,15.);
 178   fHistoTRKVtxZ->Sumw2();
 179   fOutput->Add(fHistoTRKVtxZ);
 180
 181   Int_t nBinscb=11;
 182   if(fSystem==1) nBinscb=200;
 183   if(fSystem==2) nBinscb=21;
 184   Double_t maxncn=nBinscb-0.5;
 185   fHistoNcharmed = new TH2F("hncharmed","",100,-0.5,maxMult-0.5,nBinscb,-0.5,maxncn);
 186   fHistoNcharmed->Sumw2();
 187   fOutput->Add(fHistoNcharmed);
 188   fHistoNbVsNc = new TH2F("hnbvsnc","",nBinscb,-0.5,maxncn,nBinscb,-0.5,maxncn);
 189   fHistoNbVsNc->Sumw2();
 190   fOutput->Add(fHistoNbVsNc);
 191
 192   fHistYPtPrompt[0] = new TH2F("hyptD0prompt","D0 - Prompt",40,0.,40.,20,-2.,2.);
 193   fHistYPtPrompt[1] = new TH2F("hyptDplusprompt","Dplus - Prompt",40,0.,40.,20,-2.,2.);
 194   fHistYPtPrompt[2] = new TH2F("hyptDstarprompt","Dstar - Prompt",40,0.,40.,20,-2.,2.);
 195   fHistYPtPrompt[3] = new TH2F("hyptDsprompt","Ds - Prompt",40,0.,40.,20,-2.,2.);
 196   fHistYPtPrompt[4] = new TH2F("hyptLcprompt","Lc - Prompt",40,0.,40.,20,-2.,2.);
 197
 198   fHistBYPtAllDecay[0] = new TH2F("hyptB0AllDecay","B0 - All",40,0.,40.,40,-2.,2.);
 199   fHistBYPtAllDecay[1] = new TH2F("hyptBplusAllDecay","Bplus - All",40,0.,40.,40,-2.,2.);
 200   fHistBYPtAllDecay[2] = new TH2F("hyptBstarAllDecay","Bstar - All",40,0.,40.,40,-2.,2.);
 201   fHistBYPtAllDecay[3] = new TH2F("hyptBsAllDecay","Bs - All",40,0.,40.,40,-2.,2.);
 202   fHistBYPtAllDecay[4] = new TH2F("hyptLbAllDecay","LB - All",40,0.,40.,40,-2.,2.);
 203
 204   fHistYPtAllDecay[0] = new TH2F("hyptD0AllDecay","D0 - All",40,0.,40.,40,-2.,2.);
 205   fHistYPtAllDecay[1] = new TH2F("hyptDplusAllDecay","Dplus - All",40,0.,40.,40,-2.,2.);
 206   fHistYPtAllDecay[2] = new TH2F("hyptDstarAllDecay","Dstar - All",40,0.,40.,40,-2.,2.);
 207   fHistYPtAllDecay[3] = new TH2F("hyptDsAllDecay","Ds - All",40,0.,40.,40,-2.,2.);
 208   fHistYPtAllDecay[4] = new TH2F("hyptLcAllDecay","Lc - All",40,0.,40.,40,-2.,2.);
 209
 210   fHistYPtPromptAllDecay[0] = new TH2F("hyptD0promptAllDecay","D0 - Prompt",40,0.,40.,40,-2.,2.);
 211   fHistYPtPromptAllDecay[1] = new TH2F("hyptDpluspromptAllDecay","Dplus - Prompt",40,0.,40.,40,-2.,2.);
 212   fHistYPtPromptAllDecay[2] = new TH2F("hyptDstarpromptAllDecay","Dstar - Prompt",40,0.,40.,40,-2.,2.);
 213   fHistYPtPromptAllDecay[3] = new TH2F("hyptDspromptAllDecay","Ds - Prompt",40,0.,40.,40,-2.,2.);
 214   fHistYPtPromptAllDecay[4] = new TH2F("hyptLcpromptAllDecay","Lc - Prompt",40,0.,40.,40,-2.,2.);
 215
 216   fHistYPtFeeddownAllDecay[0] = new TH2F("hyptD0feeddownAllDecay","D0 - FromB",40,0.,40.,40,-2.,2.);
 217   fHistYPtFeeddownAllDecay[1] = new TH2F("hyptDplusfeeddownAllDecay","Dplus - FromB",40,0.,40.,40,-2.,2.);
 218   fHistYPtFeeddownAllDecay[2] = new TH2F("hyptDstarfeeddownAllDecay","Dstar - FromB",40,0.,40.,40,-2.,2.);
 219   fHistYPtFeeddownAllDecay[3] = new TH2F("hyptDsfeeddownAllDecay","Ds - FromB",40,0.,40.,40,-2.,2.);
 220   fHistYPtFeeddownAllDecay[4] = new TH2F("hyptLcfeeddownAllDecay","Lc - FromB",40,0.,40.,40,-2.,2.);
 221
 222
 223   fHistYPtFeeddown[0] = new TH2F("hyptD0feeddown","D0 - Feeddown",40,0.,40.,20,-2.,2.);
 224   fHistYPtFeeddown[1] = new TH2F("hyptDplusfeeddown","Dplus - Feeddown",40,0.,40.,20,-2.,2.);
 225   fHistYPtFeeddown[2] = new TH2F("hyptDstarfeedown","Dstar - Feeddown",40,0.,40.,20,-2.,2.);
 226   fHistYPtFeeddown[3] = new TH2F("hyptDsfeedown","Ds - Feeddown",40,0.,40.,20,-2.,2.);
 227   fHistYPtFeeddown[4] = new TH2F("hyptLcfeedown","Lc - Feeddown",40,0.,40.,20,-2.,2.);
 228
 229   for(Int_t ih=0; ih<5; ih++){
 230     fHistBYPtAllDecay[ih]->Sumw2();
 231     fHistBYPtAllDecay[ih]->SetMinimum(0);
 232     fOutput->Add(fHistBYPtAllDecay[ih]);
 233     fHistYPtAllDecay[ih]->Sumw2();
 234     fHistYPtAllDecay[ih]->SetMinimum(0);
 235     fOutput->Add(fHistYPtAllDecay[ih]);
 236     fHistYPtPromptAllDecay[ih]->Sumw2();
 237     fHistYPtPromptAllDecay[ih]->SetMinimum(0);
 238     fOutput->Add(fHistYPtPromptAllDecay[ih]);
 239     fHistYPtFeeddownAllDecay[ih]->Sumw2();
 240     fHistYPtFeeddownAllDecay[ih]->SetMinimum(0);
 241     fOutput->Add(fHistYPtFeeddownAllDecay[ih]);
 242     fHistYPtPrompt[ih]->Sumw2();
 243     fHistYPtPrompt[ih]->SetMinimum(0);
 244     fOutput->Add(fHistYPtPrompt[ih]);
 245     fHistYPtFeeddown[ih]->Sumw2();
 246     fHistYPtFeeddown[ih]->SetMinimum(0);
 247     fOutput->Add(fHistYPtFeeddown[ih]);
 248   }
 249
 250   fHistYPtD0byDecChannel[0] = new TH2F("hyptD02","D0 - 2prong",40,0.,40.,20,-2.,2.);
 251   fHistYPtD0byDecChannel[1] = new TH2F("hyptD04","D0 - 4prong",40,0.,40.,20,-2.,2.);
 252   fHistYPtDplusbyDecChannel[0] = new TH2F("hyptDplusnonreson","Dplus - non reson",40,0.,40.,20,-2.,2.);
 253   fHistYPtDplusbyDecChannel[1] = new TH2F("hyptDplusreson","Dplus - reson via K0*",40,0.,40.,20,-2.,2.);
 254   fHistYPtDsbyDecChannel[0] = new TH2F("hyptDsphi","Ds - vis Phi",40,0.,40.,20,-2.,2.);
 255   fHistYPtDsbyDecChannel[1] = new TH2F("hyptDsk0st","Ds - via k0*",40,0.,40.,20,-2.,2.);
 256
 257   for(Int_t ih=0; ih<2; ih++){
 258
 259     fHistYPtD0byDecChannel[ih]->Sumw2();
 260     fHistYPtD0byDecChannel[ih]->SetMinimum(0);
 261     fOutput->Add(fHistYPtD0byDecChannel[ih]);
 262     fHistYPtDplusbyDecChannel[ih]->Sumw2();
 263     fHistYPtDplusbyDecChannel[ih]->SetMinimum(0);
 264     fOutput->Add(fHistYPtDplusbyDecChannel[ih]);
 265     fHistYPtDsbyDecChannel[ih]->Sumw2();
 266     fHistYPtDsbyDecChannel[ih]->SetMinimum(0);
 267     fOutput->Add(fHistYPtDsbyDecChannel[ih]);
 268   }
 269
 270   fHistOriginPrompt=new TH1F("hOriginPrompt","",100,0.,0.5);
 271   fHistOriginPrompt->Sumw2();
 272   fHistOriginPrompt->SetMinimum(0);
 273   fOutput->Add(fHistOriginPrompt);
 274   fHistOriginFeeddown=new TH1F("hOriginFeeddown","",100,0.,0.5);
 275   fHistOriginFeeddown->Sumw2();
 276   fHistOriginFeeddown->SetMinimum(0);
 277   fOutput->Add(fHistOriginFeeddown);
 278   fHistMotherID=new TH1F("hMotherID","",1000,-1.5,998.5);
 279   fHistMotherID->SetMinimum(0);
 280   fOutput->Add(fHistMotherID);
 281   fHistDSpecies=new TH1F("hDSpecies","",10,-0.5,9.5);
 282   fHistDSpecies->GetXaxis()->SetBinLabel(1,"D0");
 283   fHistDSpecies->GetXaxis()->SetBinLabel(2,"D0bar");
 284   fHistDSpecies->GetXaxis()->SetBinLabel(3,"D+");
 285   fHistDSpecies->GetXaxis()->SetBinLabel(4,"D-");
 286   fHistDSpecies->GetXaxis()->SetBinLabel(5,"D*+");
 287   fHistDSpecies->GetXaxis()->SetBinLabel(6,"D*-");
 288   fHistDSpecies->GetXaxis()->SetBinLabel(7,"Ds+");
 289   fHistDSpecies->GetXaxis()->SetBinLabel(8,"Ds-");
 290   fHistDSpecies->GetXaxis()->SetBinLabel(9,"Lc+");
 291   fHistDSpecies->GetXaxis()->SetBinLabel(10,"Lc-");
 292   fHistDSpecies->SetMinimum(0);
 293   fOutput->Add(fHistDSpecies);
 294   fHistBSpecies=new TH1F("hBSpecies","",10,-0.5,9.5);
 295   fHistBSpecies->GetXaxis()->SetBinLabel(1,"B0");
 296   fHistBSpecies->GetXaxis()->SetBinLabel(2,"B0bar");
 297   fHistBSpecies->GetXaxis()->SetBinLabel(3,"B+");
 298   fHistBSpecies->GetXaxis()->SetBinLabel(4,"B-");
 299   fHistBSpecies->GetXaxis()->SetBinLabel(5,"B*+");
 300   fHistBSpecies->GetXaxis()->SetBinLabel(6,"B*-");
 301   fHistBSpecies->GetXaxis()->SetBinLabel(7,"Bs+");
 302   fHistBSpecies->GetXaxis()->SetBinLabel(8,"Bs-");
 303   fHistBSpecies->GetXaxis()->SetBinLabel(9,"Lb+");
 304   fHistBSpecies->GetXaxis()->SetBinLabel(10,"Lb-");
 305   fHistBSpecies->SetMinimum(0);
 306   fOutput->Add(fHistBSpecies);
 307
 308   fHistNcollHFtype=new TH2F("hNcollHFtype","",5,-1.5,3.5,30,-0.5,29.5);
 309   fOutput->Add(fHistNcollHFtype);
 310
 311   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};
 312   const Int_t nBinsPhi=40;
 313   Double_t binsphi[nBinsPhi+1];
 314   for(Int_t ib=0; ib<=nBinsPhi; ib++) binsphi[ib]=ib*TMath::Pi()/20.;
 315   const Int_t nBinsPt=24;
 316   Double_t binspt[nBinsPt+1]={0.,0.10,0.15,0.2,0.25,
 317                               0.3,0.4,0.5,0.6,0.7,
 318                               0.8,0.9,1.,1.25,1.5,
 319                               1.75,2.,2.5,3.,4.,
 320                               5.,7.5,10.,15.,20.};
 321
 322    fHistEtaPhiPtGenEle=new TH3F("hEtaPhiPtGenEle","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
 323   fOutput->Add(fHistEtaPhiPtGenEle);
 324   fHistEtaPhiPtGenPi=new TH3F("hEtaPhiPtGenPi","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
 325   fOutput->Add(fHistEtaPhiPtGenPi);
 326   fHistEtaPhiPtGenK=new TH3F("hEtaPhiPtGenK","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
 327   fOutput->Add(fHistEtaPhiPtGenK);
 328   fHistEtaPhiPtGenPro=new TH3F("hEtaPhiPtGenPro","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
 329   fOutput->Add(fHistEtaPhiPtGenPro);
 330
 331
 332   fHistEtaPhiPtRecEle=new TH3F("hEtaPhiPtRecEle","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
 333   fOutput->Add(fHistEtaPhiPtRecEle);
 334   fHistEtaPhiPtRecPi=new TH3F("hEtaPhiPtRecPi","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
 335   fOutput->Add(fHistEtaPhiPtRecPi);
 336   fHistEtaPhiPtRecK=new TH3F("hEtaPhiPtRecK","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
 337   fOutput->Add(fHistEtaPhiPtRecK);
 338   fHistEtaPhiPtRecPro=new TH3F("hEtaPhiPtRecPro","",10,binseta,nBinsPhi,binsphi,nBinsPt,binspt);
 339   fOutput->Add(fHistEtaPhiPtRecPro);
 340
 341
 342
 343   PostData(1,fOutput);
 344
 345 }
 346 //______________________________________________________________________________
 347 void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *)
 348 {
 349   //
 350
 351   AliESDEvent *esd = (AliESDEvent*) (InputEvent());
 352
 353
 354   if(!esd) {
 355     printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
 356     return;
 357   }
 358
 359   fHistoNEvents->Fill(0);
 360
 361   if(!fESDtrackCuts){
 362     Int_t year=2011;
 363     if(esd->GetRunNumber()<=139517) year=2010;
 364     if(year==2010) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kFALSE);
 365     else fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
 366     fESDtrackCuts->SetMaxDCAToVertexXY(2.4);
 367     fESDtrackCuts->SetMaxDCAToVertexZ(3.2);
 368     fESDtrackCuts->SetDCAToVertex2D(kTRUE);
 369     fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 370                                             AliESDtrackCuts::kAny);
 371   }
 372
 373   Int_t nTracks=esd->GetNumberOfTracks();
 374   fHistoTracks->Fill(nTracks);
 375   Int_t nSelTracks=0;
 376   for(Int_t it=0; it<nTracks; it++){
 377     AliESDtrack* tr=esd->GetTrack(it);
 378     UInt_t status=tr->GetStatus();
 379     if(!(status&AliESDtrack::kITSrefit)) continue;
 380     if(!(status&AliESDtrack::kTPCin)) continue;
 381     nSelTracks++;
 382   }
 383   fHistoSelTracks->Fill(nSelTracks);
 384
 385   const AliMultiplicity* mult=esd->GetMultiplicity();
 386   Int_t nTracklets=mult->GetNumberOfTracklets();
 387   Int_t nTracklets1=0;
 388   for(Int_t it=0; it<nTracklets; it++){
 389     Double_t eta=TMath::Abs(mult->GetEta(it));
 390     if(eta<1) nTracklets1++;
 391   }
 392   fHistoTracklets->Fill(nTracklets);
 393   fHistoTrackletsEta1->Fill(nTracklets1);
 394
 395   const AliESDVertex *spdv=esd->GetVertex();
 396   if(spdv && spdv->IsFromVertexer3D()){
 397     fHistoSPD3DVtxX->Fill(spdv->GetXv());
 398     fHistoSPD3DVtxY->Fill(spdv->GetYv());
 399     fHistoSPD3DVtxZ->Fill(spdv->GetZv());
 400   }
 401   if(spdv && spdv->IsFromVertexerZ()){
 402     fHistoSPDZVtxX->Fill(spdv->GetXv());
 403     fHistoSPDZVtxY->Fill(spdv->GetYv());
 404     fHistoSPDZVtxZ->Fill(spdv->GetZv());
 405   }
 406   const AliESDVertex *trkv=esd->GetPrimaryVertex();
 407   if(trkv && trkv->GetNContributors()>1){
 408     fHistoTRKVtxX->Fill(trkv->GetXv());
 409     fHistoTRKVtxY->Fill(trkv->GetYv());
 410     fHistoTRKVtxZ->Fill(trkv->GetZv());
 411   }
 412
 413   AliStack* stack=0;
 414   if(fReadMC){
 415     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 416     if (!eventHandler) {
 417       Printf("ERROR: Could not retrieve MC event handler");
 418       return;
 419     }
 420     AliMCEvent* mcEvent = eventHandler->MCEvent();
 421     if (!mcEvent) {
 422       Printf("ERROR: Could not retrieve MC event");
 423       return;
 424     }
 425     stack = mcEvent->Stack();
 426     if (!stack) {
 427       Printf("ERROR: stack not available");
 428       return;
 429     }
 430     const AliVVertex* mcVert=mcEvent->GetPrimaryVertex();
 431     if(!mcVert){
 432       Printf("ERROR: generated vertex not available");
 433       return;
 434     }
 435     if(TMath::Abs(mcVert->GetZ())>10) return;
 436
 437     //    const AliHeader* h=(AliHeader*)mcEvent->GetHeader();
 438     //    cout<<h<<endl;
 439     TString genname=mcEvent->GenEventHeader()->ClassName();
 440     Int_t nColl=-1;
 441     Int_t typeHF=-1;
 442     TList* lgen=0x0;
 443     if(genname.Contains("CocktailEventHeader")){
 444       AliGenCocktailEventHeader *cockhead=(AliGenCocktailEventHeader*)mcEvent->GenEventHeader();
 445       lgen=cockhead->GetHeaders();
 446       for(Int_t ig=0; ig<lgen->GetEntries(); ig++){
 447         AliGenerator* gen=(AliGenerator*)lgen->At(ig);
 448         TString title=gen->GetName();
 449         if(title.Contains("bchadr")) typeHF=1;
 450         else if(title.Contains("chadr")) typeHF=0;
 451         else if(title.Contains("bele")) typeHF=3;
 452         else if(title.Contains("cele")) typeHF=2;
 453       }
 454       nColl=lgen->GetEntries();
 455       fHistNcollHFtype->Fill(typeHF,nColl);
 456     }
 457     Int_t nParticles=stack->GetNtrack();
 458     Double_t dNchdy = 0.;
 459     Int_t nb = 0, nc=0;
 460     Int_t nCharmed=0;
 461     Int_t nPhysPrim=0;
 462     for (Int_t i=0;i<nParticles;i++){
 463       TParticle* part = (TParticle*)stack->Particle(i);
 464       Int_t absPdg=TMath::Abs(part->GetPdgCode());
 465       Int_t pdg=part->GetPdgCode();
 466       if(absPdg==4) nc++;
 467       if(absPdg==5) nb++;
 468       if(stack->IsPhysicalPrimary(i)){
 469         Double_t eta=part->Eta();
 470         fHistoEtaPhysPrim->Fill(eta);
 471         if(absPdg==11) fHistEtaPhiPtGenEle->Fill(eta,part->Phi(),part->Pt());
 472         else if(absPdg==211) fHistEtaPhiPtGenPi->Fill(eta,part->Phi(),part->Pt());
 473         else if(absPdg==321) fHistEtaPhiPtGenK->Fill(eta,part->Phi(),part->Pt());
 474         else if(absPdg==2212) fHistEtaPhiPtGenPro->Fill(eta,part->Phi(),part->Pt());
 475
 476         if(TMath::Abs(eta)<0.5){
 477           dNchdy+=0.6666;   // 2/3 for the ratio charged/all
 478           nPhysPrim++;
 479         }
 480         if(TMath::Abs(eta)<0.9){
 481           fHistoPtPhysPrim->Fill(part->Pt());
 482         }
 483       }
 484       Float_t rapid=-999.;
 485       if (part->Energy() != TMath::Abs(part->Pz())){
 486         rapid=0.5*TMath::Log((part->Energy()+part->Pz())/(part->Energy()-part->Pz()));
 487       }
 488
 489
 490       Int_t iPart=-1;
 491       Int_t iType=0;
 492       Int_t iSpecies=-1;
 493       if(absPdg==421){
 494         iSpecies=0;
 495         iType=CheckD0Decay(i,stack);
 496         if(iType>=0) iPart=0;
 497       }
 498       else if(absPdg==411){
 499         iSpecies=1;
 500         iType=CheckDplusDecay(i,stack);
 501         if(iType>=0) iPart=1;
 502       }
 503       else if(absPdg==413){
 504         iSpecies=2;
 505         iType=CheckDstarDecay(i,stack);
 506         if(iType>=0) iPart=2;
 507       }
 508       else if(absPdg==431){
 509         iSpecies=3;
 510         iType=CheckDsDecay(i,stack);
 511         if(iType==0 || iType==1) iPart=3;
 512       }
 513       else if(absPdg==4122){
 514         iSpecies=4;
 515         iType=CheckLcDecay(i,stack);
 516         if(iType>=0) iPart=4;
 517       }
 518       if(iSpecies>=0) fHistYPtAllDecay[iSpecies]->Fill(part->Pt(),rapid);
 519
 520       // check beauty mesons
 521       if(absPdg==511) fHistBYPtAllDecay[0]->Fill(part->Pt(),rapid);
 522       else if(absPdg==521) fHistBYPtAllDecay[1]->Fill(part->Pt(),rapid);
 523       else if(absPdg==513) fHistBYPtAllDecay[2]->Fill(part->Pt(),rapid);
 524       else if(absPdg==531) fHistBYPtAllDecay[3]->Fill(part->Pt(),rapid);
 525       else if(absPdg==5122) fHistBYPtAllDecay[4]->Fill(part->Pt(),rapid);
 526
 527       if(pdg==511) fHistBSpecies->Fill(0);
 528       else if(pdg==-511) fHistBSpecies->Fill(1);
 529       else if(pdg==521) fHistBSpecies->Fill(2);
 530       else if(pdg==-521) fHistBSpecies->Fill(3);
 531       else if(pdg==513) fHistBSpecies->Fill(4);
 532       else if(pdg==-513) fHistBSpecies->Fill(5);
 533       else if(pdg==531) fHistBSpecies->Fill(6);
 534       else if(pdg==-531) fHistBSpecies->Fill(7);
 535       else if(pdg==5122) fHistBSpecies->Fill(8);
 536       else if(pdg==-5122) fHistBSpecies->Fill(9);
 537
 538      if(iSpecies<0) continue; // not a charm meson
 539
 540       if(pdg==421) fHistDSpecies->Fill(0);
 541       else if(pdg==-421) fHistDSpecies->Fill(1);
 542       else if(pdg==411) fHistDSpecies->Fill(2);
 543       else if(pdg==-411) fHistDSpecies->Fill(3);
 544       else if(pdg==413) fHistDSpecies->Fill(4);
 545       else if(pdg==-413) fHistDSpecies->Fill(5);
 546       else if(pdg==431) fHistDSpecies->Fill(6);
 547       else if(pdg==-431) fHistDSpecies->Fill(7);
 548       else if(pdg==4122) fHistDSpecies->Fill(8);
 549       else if(pdg==-4122) fHistDSpecies->Fill(9);
 550
 551       Double_t distx=part->Vx()-mcVert->GetX();
 552       Double_t disty=part->Vy()-mcVert->GetY();
 553       Double_t distz=part->Vz()-mcVert->GetZ();
 554       Double_t distToVert=TMath::Sqrt(distx*distx+disty*disty+distz*distz);
 555       fHistMotherID->Fill(part->GetFirstMother());
 556       TParticle* runningpart=part;
 557       Int_t iFromB=-1;
 558       Int_t pdgmoth=-1;
 559       if(fSearchUpToQuark){
 560         while(1){
 561           Int_t labmoth=runningpart->GetFirstMother();
 562           if(labmoth==-1) break;
 563           TParticle *mot=(TParticle*)stack->Particle(labmoth);
 564           pdgmoth=TMath::Abs(mot->GetPdgCode());
 565           if(pdgmoth==5){
 566             iFromB=1;
 567             break;
 568           }else if(pdgmoth==4){
 569             iFromB=0;
 570             break;
 571           }
 572           runningpart=mot;
 573         }
 574       }else{
 575         iFromB=0;
 576         while(1){
 577           Int_t labmoth=runningpart->GetFirstMother();
 578           if(labmoth==-1) break;
 579           TParticle *mot=(TParticle*)stack->Particle(labmoth);
 580           pdgmoth=TMath::Abs(mot->GetPdgCode());
 581           if(pdgmoth>=500 && pdgmoth<=599){
 582             iFromB=1;
 583             break;
 584           }
 585           if(pdgmoth>=5000 && pdgmoth<=5999){
 586             iFromB=1;
 587             break;
 588           }
 589           runningpart=mot;
 590         }
 591       }
 592
 593       if(iFromB==0){
 594         fHistYPtPromptAllDecay[iSpecies]->Fill(part->Pt(),rapid);
 595         fHistOriginPrompt->Fill(distToVert);
 596       }
 597       else if(iFromB==1){
 598         fHistYPtFeeddownAllDecay[iSpecies]->Fill(part->Pt(),rapid);
 599         fHistOriginFeeddown->Fill(distToVert);
 600       }
 601
 602       if(iPart<0) continue;
 603       if(iType<0) continue;
 604       nCharmed++;
 605       if(iPart==0 && iType<=1){
 606         fHistYPtD0byDecChannel[iType]->Fill(part->Pt(),rapid);
 607       }else if(iPart==1 && iType<=1){
 608         fHistYPtDplusbyDecChannel[iType]->Fill(part->Pt(),rapid);
 609       }else if(iPart==3 && iType<=1){
 610         fHistYPtDsbyDecChannel[iType]->Fill(part->Pt(),rapid);
 611       }
 612
 613       if(iFromB==0 && iPart>=0 && iPart<5) fHistYPtPrompt[iPart]->Fill(part->Pt(),rapid);
 614       else if(iFromB==1 && iPart>=0 && iPart<5) fHistYPtFeeddown[iPart]->Fill(part->Pt(),rapid);
 615     }
 616
 617     for(Int_t i=0; i<nTracks; i++){
 618       AliESDtrack* track=esd->GetTrack(i);
 619       if(fESDtrackCuts->AcceptTrack(track)){
 620         Int_t label=TMath::Abs(track->GetLabel());
 621         if(stack->IsPhysicalPrimary(label)){
 622           TParticle* part = (TParticle*)stack->Particle(label);
 623           Int_t absPdg=TMath::Abs(part->GetPdgCode());
 624           Double_t eta=part->Eta();
 625           if(absPdg==11) fHistEtaPhiPtRecEle->Fill(eta,part->Phi(),part->Pt());
 626           else if(absPdg==211) fHistEtaPhiPtRecPi->Fill(eta,part->Phi(),part->Pt());
 627           else if(absPdg==321) fHistEtaPhiPtRecK->Fill(eta,part->Phi(),part->Pt());
 628           else if(absPdg==2212) fHistEtaPhiPtRecPro->Fill(eta,part->Phi(),part->Pt());
 629         }
 630       }
 631     }
 632     fHistoNcharmed->Fill(dNchdy,nCharmed);
 633     fHistoNbVsNc->Fill(nc,nb);
 634     fHistoPhysPrim->Fill(nPhysPrim);
 635   }
 636
 637   PostData(1,fOutput);
 638
 639 }
 640 //______________________________________________________________________________
 641 void AliAnalysisTaskCheckHFMCProd::Terminate(Option_t */*option*/)
 642 {
 643   // Terminate analysis
 644   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 645   if (!fOutput) {
 646     printf("ERROR: fOutput not available\n");
 647     return;
 648   }
 649
 650   return;
 651 }
 652
 653
 654
 655
 656 //______________________________________________________________________________
 657 Int_t AliAnalysisTaskCheckHFMCProd::CheckD0Decay(Int_t labD0, AliStack* stack) const{
 658
 659   if(labD0<0) return -1;
 660   TParticle* dp = (TParticle*)stack->Particle(labD0);
 661   Int_t pdgdp=dp->GetPdgCode();
 662   Int_t nDau=dp->GetNDaughters();
 663
 664   if(nDau==2){
 665     Int_t nKaons=0;
 666     Int_t nPions=0;
 667     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 668       if(iDau<0) return -1;
 669       TParticle* dau=(TParticle*)stack->Particle(iDau);
 670       Int_t pdgdau=dau->GetPdgCode();
 671       if(TMath::Abs(pdgdau)==321){
 672         if(pdgdp>0 && pdgdau>0) return -1;
 673         if(pdgdp<0 && pdgdau<0) return -1;
 674         nKaons++;
 675       }else if(TMath::Abs(pdgdau)==211){
 676         if(pdgdp<0 && pdgdau>0) return -1;
 677         if(pdgdp>0 && pdgdau<0) return -1;
 678         nPions++;
 679       }else{
 680         return -1;
 681       }
 682     }
 683     if(nPions!=1) return -1;
 684     if(nKaons!=1) return -1;
 685     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 686       if(iDau<0) return -1;
 687     }
 688     return 0;
 689   }
 690
 691   if(nDau==3 || nDau==4){
 692     Int_t nKaons=0;
 693     Int_t nPions=0;
 694     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 695       if(iDau<0) return -1;
 696       TParticle* dau=(TParticle*)stack->Particle(iDau);
 697       Int_t pdgdau=dau->GetPdgCode();
 698       if(TMath::Abs(pdgdau)==321){
 699         if(pdgdp>0 && pdgdau>0) return -1;
 700         if(pdgdp<0 && pdgdau<0) return -1;
 701         nKaons++;
 702       }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
 703         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 704           if(resDau<0) return -1;
 705           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 706           Int_t pdgresdau=resdau->GetPdgCode();
 707           if(TMath::Abs(pdgresdau)==321){
 708             if(pdgdp>0 && pdgresdau>0) return -1;
 709             if(pdgdp<0 && pdgresdau<0) return -1;
 710             nKaons++;
 711           }
 712           if(TMath::Abs(pdgresdau)==211){
 713             nPions++;
 714           }
 715         }
 716       }else if(TMath::Abs(pdgdau)==211){
 717           nPions++;
 718       }else{
 719         return -1;
 720       }
 721     }
 722     if(nPions!=3) return -1;
 723     if(nKaons!=1) return -1;
 724     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 725       if(iDau<0) return -1;
 726     }
 727     return 1;
 728   }
 729
 730   return -1;
 731 }
 732
 733
 734 //______________________________________________________________________________
 735 Int_t AliAnalysisTaskCheckHFMCProd::CheckDplusDecay(Int_t labDplus, AliStack* stack) const{
 736
 737   if(labDplus<0) return -1;
 738   TParticle* dp = (TParticle*)stack->Particle(labDplus);
 739   Int_t pdgdp=dp->GetPdgCode();
 740   Int_t nDau=dp->GetNDaughters();
 741
 742   if(nDau==3){
 743     Int_t nKaons=0;
 744     Int_t nPions=0;
 745     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 746       if(iDau<0) return -1;
 747       TParticle* dau=(TParticle*)stack->Particle(iDau);
 748       Int_t pdgdau=dau->GetPdgCode();
 749       if(TMath::Abs(pdgdau)==321){
 750         if(pdgdp>0 && pdgdau>0) return -1;
 751         if(pdgdp<0 && pdgdau<0) return -1;
 752         nKaons++;
 753       }else if(TMath::Abs(pdgdau)==211){
 754         if(pdgdp<0 && pdgdau>0) return -1;
 755         if(pdgdp>0 && pdgdau<0) return -1;
 756         nPions++;
 757       }else{
 758         return -1;
 759       }
 760     }
 761     if(nPions!=2) return -1;
 762     if(nKaons!=1) return -1;
 763     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 764       if(iDau<0) return -1;
 765     }
 766     return 0;
 767   }
 768
 769   if(nDau==2){
 770     Int_t nKaons=0;
 771     Int_t nPions=0;
 772     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 773       if(iDau<0) return -1;
 774       TParticle* dau=(TParticle*)stack->Particle(iDau);
 775       Int_t pdgdau=dau->GetPdgCode();
 776       if(TMath::Abs(pdgdau)==313){
 777         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 778           if(resDau<0) return -1;
 779           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 780           Int_t pdgresdau=resdau->GetPdgCode();
 781           if(TMath::Abs(pdgresdau)==321){
 782             if(pdgdp>0 && pdgresdau>0) return -1;
 783             if(pdgdp<0 && pdgresdau<0) return -1;
 784             nKaons++;
 785           }
 786           if(TMath::Abs(pdgresdau)==211){
 787             if(pdgdp<0 && pdgresdau>0) return -1;
 788             if(pdgdp>0 && pdgresdau<0) return -1;
 789             nPions++;
 790           }
 791         }
 792       }else if(TMath::Abs(pdgdau)==211){
 793           if(pdgdp<0 && pdgdau>0) return -1;
 794           if(pdgdp>0 && pdgdau<0) return -1;
 795           nPions++;
 796       }else{
 797         return -1;
 798       }
 799     }
 800     if(nPions!=2) return -1;
 801     if(nKaons!=1) return -1;
 802     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 803       if(iDau<0) return -1;
 804     }
 805     return 1;
 806   }
 807   return -1;
 808 }
 809
 810 //______________________________________________________________________________
 811 Int_t AliAnalysisTaskCheckHFMCProd::CheckDsDecay(Int_t labDs, AliStack* stack) const{
 812   // Ds decay
 813   if(labDs<0) return -1;
 814   TParticle* dp = (TParticle*)stack->Particle(labDs);
 815   Int_t pdgdp=dp->GetPdgCode();
 816   Int_t nDau=dp->GetNDaughters();
 817
 818   if(nDau==3){
 819     Int_t nKaons=0;
 820     Int_t nPions=0;
 821     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 822       if(iDau<0) return -1;
 823       TParticle* dau=(TParticle*)stack->Particle(iDau);
 824       Int_t pdgdau=dau->GetPdgCode();
 825       if(TMath::Abs(pdgdau)==321){
 826         nKaons++;
 827       }else if(TMath::Abs(pdgdau)==211){
 828         if(pdgdp<0 && pdgdau>0) return -1;
 829         if(pdgdp>0 && pdgdau<0) return -1;
 830         nPions++;
 831       }else{
 832         return -1;
 833       }
 834     }
 835     if(nPions!=1) return -1;
 836     if(nKaons!=2) return -1;
 837     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 838       if(iDau<0) return -1;
 839     }
 840     return 2;
 841   }
 842
 843   if(nDau==2){
 844     Int_t nKaons=0;
 845     Int_t nPions=0;
 846     Bool_t isPhi=kFALSE;
 847     Bool_t isk0st=kFALSE;
 848     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 849       if(iDau<0) return -1;
 850       TParticle* dau=(TParticle*)stack->Particle(iDau);
 851       Int_t pdgdau=dau->GetPdgCode();
 852       if(TMath::Abs(pdgdau)==313){
 853         isk0st=kTRUE;
 854         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 855           if(resDau<0) return -1;
 856           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 857           Int_t pdgresdau=resdau->GetPdgCode();
 858           if(TMath::Abs(pdgresdau)==321){
 859             nKaons++;
 860           }
 861           if(TMath::Abs(pdgresdau)==211){
 862             if(pdgdp<0 && pdgresdau>0) return -1;
 863             if(pdgdp>0 && pdgresdau<0) return -1;
 864             nPions++;
 865           }
 866         }
 867       }else if(TMath::Abs(pdgdau)==333){
 868         isPhi=kTRUE;
 869         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 870           if(resDau<0) return -1;
 871           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 872           if(!resdau) return -1;
 873           Int_t pdgresdau=resdau->GetPdgCode();
 874           if(TMath::Abs(pdgresdau)==321){
 875             nKaons++;
 876           }else{
 877             return -1;
 878           }
 879         }
 880       }else if(TMath::Abs(pdgdau)==211){
 881         if(pdgdp<0 && pdgdau>0) return -1;
 882         if(pdgdp>0 && pdgdau<0) return -1;
 883         nPions++;
 884       }else if(TMath::Abs(pdgdau)==321){
 885         nKaons++;
 886       }else{
 887         return -1;
 888       }
 889     }
 890     if(nPions!=1) return -1;
 891     if(nKaons!=2) return -1;
 892     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 893       if(iDau<0) return -1;
 894     }
 895     if(isk0st) return 1;
 896     else if(isPhi) return 0;
 897     else return 3;
 898   }
 899   return -1;
 900 }
 901
 902 //______________________________________________________________________________
 903 Int_t AliAnalysisTaskCheckHFMCProd::CheckDstarDecay(Int_t labDstar, AliStack* stack) const{
 904
 905   if(labDstar<0) return -1;
 906   TParticle* dp = (TParticle*)stack->Particle(labDstar);
 907   Int_t pdgdp=dp->GetPdgCode();
 908   Int_t nDau=dp->GetNDaughters();
 909   if(nDau!=2) return -1;
 910
 911   Int_t nKaons=0;
 912   Int_t nPions=0;
 913   for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 914     if(iDau<0) return -1;
 915     TParticle* dau=(TParticle*)stack->Particle(iDau);
 916     Int_t pdgdau=dau->GetPdgCode();
 917     if(TMath::Abs(pdgdau)==421){
 918       for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 919         if(resDau<0) return -1;
 920         TParticle* resdau=(TParticle*)stack->Particle(resDau);
 921         Int_t pdgresdau=resdau->GetPdgCode();
 922         if(TMath::Abs(pdgresdau)==321){
 923           if(pdgdp>0 && pdgresdau>0) return -1;
 924           if(pdgdp<0 && pdgresdau<0) return -1;
 925           nKaons++;
 926         }
 927         if(TMath::Abs(pdgresdau)==211){
 928           if(pdgdp<0 && pdgresdau>0) return -1;
 929           if(pdgdp>0 && pdgresdau<0) return -1;
 930           nPions++;
 931         }
 932       }
 933     }else if(TMath::Abs(pdgdau)==211){
 934       if(pdgdp<0 && pdgdau>0) return -1;
 935       if(pdgdp>0 && pdgdau<0) return -1;
 936       nPions++;
 937     }else{
 938       return -1;
 939     }
 940   }
 941   if(nPions!=2) return -1;
 942   if(nKaons!=1) return -1;
 943   for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 944     if(iDau<0) return -1;
 945   }
 946   return 0;
 947
 948 }
 949
 950 //______________________________________________________________________________
 951 Int_t AliAnalysisTaskCheckHFMCProd::CheckLcDecay(Int_t labLc, AliStack* stack) const{
 952   if(labLc<0) return -1;
 953   TParticle* dp = (TParticle*)stack->Particle(labLc);
 954   Int_t pdgdp=dp->GetPdgCode();
 955   Int_t nDau=dp->GetNDaughters();
 956
 957   if(nDau==3){
 958     Int_t nKaons=0;
 959     Int_t nPions=0;
 960     Int_t nProtons=0;
 961     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 962       if(iDau<0) return -1;
 963       TParticle* dau=(TParticle*)stack->Particle(iDau);
 964       Int_t pdgdau=dau->GetPdgCode();
 965       if(TMath::Abs(pdgdau)==321){
 966         if(pdgdp>0 && pdgdau>0) return -1;
 967         if(pdgdp<0 && pdgdau<0) return -1;
 968         nKaons++;
 969       }else if(TMath::Abs(pdgdau)==211){
 970         if(pdgdp<0 && pdgdau>0) return -1;
 971         if(pdgdp>0 && pdgdau<0) return -1;
 972         nPions++;
 973       }else if(TMath::Abs(pdgdau)==2212){
 974         if(pdgdp<0 && pdgdau>0) return -1;
 975         if(pdgdp>0 && pdgdau<0) return -1;
 976         nProtons++;
 977       }else{
 978         return -1;
 979       }
 980     }
 981     if(nPions!=1) return -1;
 982     if(nKaons!=1) return -1;
 983     if(nProtons!=1) return -1;
 984     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 985       if(iDau<0) return -1;
 986     }
 987     return 0;
 988   }
 989
 990   if(nDau==2){
 991     Int_t nKaons=0;
 992     Int_t nPions=0;
 993     Int_t nProtons=0;
 994     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 995       if(iDau<0) return -1;
 996       TParticle* dau=(TParticle*)stack->Particle(iDau);
 997       Int_t pdgdau=dau->GetPdgCode();
 998       if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 || TMath::Abs(pdgdau)==2224
 999          || TMath::Abs(pdgdau)==3122 || TMath::Abs(pdgdau)==311){
1000         Int_t nDauRes=dau->GetNDaughters();
1001         if(nDauRes!=2)  return -1;
1002         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
1003           if(resDau<0) return -1;
1004           TParticle* resdau=(TParticle*)stack->Particle(resDau);
1005           Int_t pdgresdau=resdau->GetPdgCode();
1006           if(TMath::Abs(pdgresdau)==321){
1007             if(pdgdp>0 && pdgresdau>0) return -1;
1008             if(pdgdp<0 && pdgresdau<0) return -1;
1009             nKaons++;
1010           }
1011           if(TMath::Abs(pdgresdau)==211){
1012             if(pdgdp<0 && pdgresdau>0) return -1;
1013             if(pdgdp>0 && pdgresdau<0) return -1;
1014             nPions++;
1015           }
1016           if(TMath::Abs(pdgresdau)==2212){
1017             if(pdgdp<0 && pdgresdau>0) return -1;
1018             if(pdgdp>0 && pdgresdau<0) return -1;
1019             nProtons++;
1020           }
1021         }
1022       }else if(TMath::Abs(pdgdau)==321){
1023         if(pdgdp>0 && pdgdau>0) return -1;
1024         if(pdgdp<0 && pdgdau<0) return -1;
1025         nKaons++;
1026       }else if(TMath::Abs(pdgdau)==211){
1027         if(pdgdp<0 && pdgdau>0) return -1;
1028         if(pdgdp>0 && pdgdau<0) return -1;
1029         nPions++;
1030       }else if(TMath::Abs(pdgdau)==2212){
1031         if(pdgdp<0 && pdgdau>0) return -1;
1032         if(pdgdp>0 && pdgdau<0) return -1;
1033         nProtons++;
1034       }else{
1035         return -1;
1036       }
1037     }
1038     if(nPions!=1) return -1;
1039     if(nKaons!=1) return -1;
1040     if(nProtons!=1) return -1;
1041     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
1042       if(iDau<0) return -1;
1043     }
1044     return 1;
1045   }
1046   return -1;
1047 }
1048