PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.cxx

   1 #include "AliAnalysisTaskSE.h"
   2 #include "AliAnalysisManager.h"
   3 #include "AliAnalysisDataContainer.h"
   4 #include "AliESDEvent.h"
   5 #include "AliStack.h"
   6 #include "AliCentrality.h"
   7 #include "AliMCEventHandler.h"
   8 #include "AliMCEvent.h"
   9 #include "AliMultiplicity.h"
  10 #include <TParticle.h>
  11 #include <TSystem.h>
  12 #include <TTree.h>
  13 #include <TNtuple.h>
  14 #include <TH1F.h>
  15 #include <TH2F.h>
  16 #include <TChain.h>
  17 #include "AliESDInputHandlerRP.h"
  18 #include "AliAnalysisTaskCheckHFMCProd.h"
  19
  20 /**************************************************************************
  21  * Copyright(c) 1998-2012, ALICE Experiment at CERN, All rights reserved. *
  22  *                                                                        *
  23  * Author: The ALICE Off-line Project.                                    *
  24  * Contributors are mentioned in the code where appropriate.              *
  25  *                                                                        *
  26  * Permission to use, copy, modify and distribute this software and its   *
  27  * documentation strictly for non-commercial purposes is hereby granted   *
  28  * without fee, provided that the above copyright notice appears in all   *
  29  * copies and that both the copyright notice and this permission notice   *
  30  * appear in the supporting documentation. The authors make no claims     *
  31  * about the suitability of this software for any purpose. It is          *
  32  * provided "as is" without express or implied warranty.                  *
  33  **************************************************************************/
  34
  35 /* $Id$ */
  36
  37 //*************************************************************************
  38 // Implementation of class AliAnalysisTaskCheckHFMCProd
  39 // AliAnalysisTask to check MC production at ESD+Kine level
  40 //
  41 //
  42 // Authors: F. Prino, prino@to.infn.it
  43 //
  44 //*************************************************************************
  45
  46 ClassImp(AliAnalysisTaskCheckHFMCProd)
  47 //______________________________________________________________________________
  48 AliAnalysisTaskCheckHFMCProd::AliAnalysisTaskCheckHFMCProd() : AliAnalysisTaskSE("HFMCChecks"),
  49   fOutput(0),
  50   fHistoNEvents(0),
  51   fHistoTracks(0),
  52   fHistoSelTracks(0),
  53   fHistoTracklets(0),
  54   fHistoSPD3DVtxX(0),
  55   fHistoSPD3DVtxY(0),
  56   fHistoSPD3DVtxZ(0),
  57   fHistoSPDZVtxX(0),
  58   fHistoSPDZVtxY(0),
  59   fHistoSPDZVtxZ(0),
  60   fHistoTRKVtxX(0),
  61   fHistoTRKVtxY(0),
  62   fHistoTRKVtxZ(0),
  63   fHistoNcharmed(0),
  64   fHistoNbVsNc(0),
  65   fPbPb(kFALSE),
  66   fReadMC(kTRUE)
  67 {
  68   //
  69   DefineInput(0, TChain::Class());
  70   DefineOutput(1, TList::Class());
  71 }
  72
  73
  74 //___________________________________________________________________________
  75 AliAnalysisTaskCheckHFMCProd::~AliAnalysisTaskCheckHFMCProd(){
  76   //
  77   if (AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return;
  78   if (fOutput) {
  79     delete fOutput;
  80     fOutput = 0;
  81   }
  82 }
  83
  84 //___________________________________________________________________________
  85 void AliAnalysisTaskCheckHFMCProd::UserCreateOutputObjects() {
  86   // create output histos
  87
  88   fOutput = new TList();
  89   fOutput->SetOwner();
  90   fOutput->SetName("OutputHistos");
  91
  92   fHistoNEvents = new TH1F("hNEvents", "Number of processed events",3,-0.5,2.5);
  93   fHistoNEvents->Sumw2();
  94   fHistoNEvents->SetMinimum(0);
  95   fOutput->Add(fHistoNEvents);
  96
  97   Double_t maxMult=100.;
  98   if(fPbPb) maxMult=10000.;
  99   fHistoTracks = new TH1F("hTracks","",100,0.,maxMult*2);
 100   fHistoTracks->Sumw2();
 101   fOutput->Add(fHistoTracks);
 102   fHistoSelTracks = new TH1F("hSelTracks","",100,0.,maxMult);
 103   fHistoSelTracks->Sumw2();
 104   fOutput->Add(fHistoSelTracks);
 105   fHistoTracklets = new TH1F("hTracklets","",100,0.,maxMult);
 106   fHistoTracklets->Sumw2();
 107   fOutput->Add(fHistoTracklets);
 108
 109   fHistoSPD3DVtxX = new TH1F("hSPD3DvX","",100,-1.,1.);
 110   fHistoSPD3DVtxX->Sumw2();
 111   fOutput->Add(fHistoSPD3DVtxX);
 112   fHistoSPD3DVtxY = new TH1F("hSPD3DvY","",100,-1.,1.);
 113   fHistoSPD3DVtxY->Sumw2();
 114   fOutput->Add(fHistoSPD3DVtxY);
 115   fHistoSPD3DVtxZ = new TH1F("hSPD3DvZ","",100,-15.,15.);
 116   fHistoSPD3DVtxZ->Sumw2();
 117   fOutput->Add(fHistoSPD3DVtxZ);
 118
 119   fHistoSPDZVtxX = new TH1F("hSPDZvX","",100,-1.,1.);
 120   fHistoSPDZVtxX->Sumw2();
 121   fOutput->Add(fHistoSPDZVtxX);
 122   fHistoSPDZVtxY = new TH1F("hSPDZvY","",100,-1.,1.);
 123   fHistoSPDZVtxY->Sumw2();
 124   fOutput->Add(fHistoSPDZVtxY);
 125   fHistoSPDZVtxZ = new TH1F("hSPDZvZ","",100,-15.,15.);
 126   fHistoSPDZVtxZ->Sumw2();
 127   fOutput->Add(fHistoSPDZVtxZ);
 128
 129
 130   fHistoTRKVtxX = new TH1F("hTRKvX","",100,-1.,1.);
 131   fHistoTRKVtxX->Sumw2();
 132   fOutput->Add(fHistoTRKVtxX);
 133   fHistoTRKVtxY = new TH1F("hTRKvY","",100,-1.,1.);
 134   fHistoTRKVtxY->Sumw2();
 135   fOutput->Add(fHistoTRKVtxY);
 136   fHistoTRKVtxZ = new TH1F("hTRKvZ","",100,-15.,15.);
 137   fHistoTRKVtxZ->Sumw2();
 138   fOutput->Add(fHistoTRKVtxZ);
 139
 140   Int_t nBinscb=11;
 141   if(fPbPb) nBinscb=200;
 142   Double_t maxncn=nBinscb-0.5;
 143   fHistoNcharmed = new TH2F("hncharmed","",100,0.,maxMult,nBinscb,-0.5,maxncn);
 144   fHistoNcharmed->Sumw2();
 145   fOutput->Add(fHistoNcharmed);
 146   fHistoNbVsNc = new TH2F("hnbvsnc","",nBinscb,-0.5,maxncn,nBinscb,-0.5,maxncn);
 147   fHistoNbVsNc->Sumw2();
 148   fOutput->Add(fHistoNbVsNc);
 149
 150   fHistYPtPrompt[0] = new TH2F("hyptd0prompt","D0 - Prompt",20,0.,20.,20,-2.,2.);
 151   fHistYPtPrompt[1] = new TH2F("hyptdplusprompt","Dplus - Prompt",20,0.,20.,20,-2.,2.);
 152   fHistYPtPrompt[2] = new TH2F("hyptdstarprompt","Dstar - Prompt",20,0.,20.,20,-2.,2.);
 153   fHistYPtPrompt[3] = new TH2F("hyptdsprompt","Ds - Prompt",20,0.,20.,20,-2.,2.);
 154   fHistYPtPrompt[4] = new TH2F("hyptlcprompt","Lc - Prompt",20,0.,20.,20,-2.,2.);
 155
 156   fHistYPtAllDecay[0] = new TH2F("hyptd0AllDecay","D0 - All",20,0.,20.,40,-2.,2.);
 157   fHistYPtAllDecay[1] = new TH2F("hyptdplusAllDecay","Dplus - All",20,0.,20.,40,-2.,2.);
 158   fHistYPtAllDecay[2] = new TH2F("hyptdstarAllDecay","Dstar - All",20,0.,20.,40,-2.,2.);
 159   fHistYPtAllDecay[3] = new TH2F("hyptdsAllDecay","Ds - All",20,0.,20.,40,-2.,2.);
 160   fHistYPtAllDecay[4] = new TH2F("hyptlcAllDecay","Lc - All",20,0.,20.,40,-2.,2.);
 161
 162   fHistYPtPromptAllDecay[0] = new TH2F("hyptd0promptAllDecay","D0 - Prompt",20,0.,20.,40,-2.,2.);
 163   fHistYPtPromptAllDecay[1] = new TH2F("hyptdpluspromptAllDecay","Dplus - Prompt",20,0.,20.,40,-2.,2.);
 164   fHistYPtPromptAllDecay[2] = new TH2F("hyptdstarpromptAllDecay","Dstar - Prompt",20,0.,20.,40,-2.,2.);
 165   fHistYPtPromptAllDecay[3] = new TH2F("hyptdspromptAllDecay","Ds - Prompt",20,0.,20.,40,-2.,2.);
 166   fHistYPtPromptAllDecay[4] = new TH2F("hyptlcpromptAllDecay","Lc - Prompt",20,0.,20.,40,-2.,2.);
 167
 168   fHistYPtFeeddownAllDecay[0] = new TH2F("hyptd0feeddownAllDecay","D0 - FromB",20,0.,20.,40,-2.,2.);
 169   fHistYPtFeeddownAllDecay[1] = new TH2F("hyptdplusfeeddownAllDecay","Dplus - FromB",20,0.,20.,40,-2.,2.);
 170   fHistYPtFeeddownAllDecay[2] = new TH2F("hyptdstarfeeddownAllDecay","Dstar - FromB",20,0.,20.,40,-2.,2.);
 171   fHistYPtFeeddownAllDecay[3] = new TH2F("hyptdsfeeddownAllDecay","Ds - FromB",20,0.,20.,40,-2.,2.);
 172   fHistYPtFeeddownAllDecay[4] = new TH2F("hyptlcfeeddownAllDecay","Lc - FromB",20,0.,20.,40,-2.,2.);
 173
 174
 175  fHistYPtFeeddown[0] = new TH2F("hyptd0feeddown","D0 - Feeddown",20,0.,20.,20,-2.,2.);
 176   fHistYPtFeeddown[1] = new TH2F("hyptdplusfeeddown","Dplus - Feeddown",20,0.,20.,20,-2.,2.);
 177   fHistYPtFeeddown[2] = new TH2F("hyptdstarfeedown","Dstar - Feeddown",20,0.,20.,20,-2.,2.);
 178   fHistYPtFeeddown[3] = new TH2F("hyptdsfeedown","Ds - Feeddown",20,0.,20.,20,-2.,2.);
 179   fHistYPtFeeddown[4] = new TH2F("hyptlcfeedown","Lc - Feeddown",20,0.,20.,20,-2.,2.);
 180
 181   for(Int_t ih=0; ih<5; ih++){
 182     fHistYPtAllDecay[ih]->Sumw2();
 183     fHistYPtAllDecay[ih]->SetMinimum(0);
 184     fOutput->Add(fHistYPtAllDecay[ih]);
 185     fHistYPtPromptAllDecay[ih]->Sumw2();
 186     fHistYPtPromptAllDecay[ih]->SetMinimum(0);
 187     fOutput->Add(fHistYPtPromptAllDecay[ih]);
 188     fHistYPtFeeddownAllDecay[ih]->Sumw2();
 189     fHistYPtFeeddownAllDecay[ih]->SetMinimum(0);
 190     fOutput->Add(fHistYPtFeeddownAllDecay[ih]);
 191     fHistYPtPrompt[ih]->Sumw2();
 192     fHistYPtPrompt[ih]->SetMinimum(0);
 193     fOutput->Add(fHistYPtPrompt[ih]);
 194     fHistYPtFeeddown[ih]->Sumw2();
 195     fHistYPtFeeddown[ih]->SetMinimum(0);
 196     fOutput->Add(fHistYPtFeeddown[ih]);
 197   }
 198
 199   fHistYPtD0byDecChannel[0] = new TH2F("hyptD02","D0 - 2prong",20,0.,20.,20,-2.,2.);
 200   fHistYPtD0byDecChannel[1] = new TH2F("hyptD04","D0 - 4prong",20,0.,20.,20,-2.,2.);
 201   fHistYPtDplusbyDecChannel[0] = new TH2F("hyptDplusnonreson","Dplus - non reson",20,0.,20.,20,-2.,2.);
 202   fHistYPtDplusbyDecChannel[1] = new TH2F("hyptDplusreson","Dplus - reson via K0*",20,0.,20.,20,-2.,2.);
 203   fHistYPtDsbyDecChannel[0] = new TH2F("hyptdsphi","Ds - vis Phi",20,0.,20.,20,-2.,2.);
 204   fHistYPtDsbyDecChannel[1] = new TH2F("hyptdsk0st","Ds - via k0*",20,0.,20.,20,-2.,2.);
 205
 206   for(Int_t ih=0; ih<2; ih++){
 207
 208     fHistYPtD0byDecChannel[ih]->Sumw2();
 209     fHistYPtD0byDecChannel[ih]->SetMinimum(0);
 210     fOutput->Add(fHistYPtD0byDecChannel[ih]);
 211     fHistYPtDplusbyDecChannel[ih]->Sumw2();
 212     fHistYPtDplusbyDecChannel[ih]->SetMinimum(0);
 213     fOutput->Add(fHistYPtDplusbyDecChannel[ih]);
 214     fHistYPtDsbyDecChannel[ih]->Sumw2();
 215     fHistYPtDsbyDecChannel[ih]->SetMinimum(0);
 216     fOutput->Add(fHistYPtDsbyDecChannel[ih]);
 217   }
 218
 219
 220   PostData(1,fOutput);
 221
 222 }
 223 //______________________________________________________________________________
 224 void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *)
 225 {
 226   //
 227
 228   AliESDEvent *esd = (AliESDEvent*) (InputEvent());
 229
 230
 231   if(!esd) {
 232     printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
 233     return;
 234   }
 235
 236   fHistoNEvents->Fill(0);
 237
 238   Int_t nTracks=esd->GetNumberOfTracks();
 239   fHistoTracks->Fill(nTracks);
 240   Int_t nSelTracks=0;
 241   for(Int_t it=0; it<nTracks; it++){
 242     AliESDtrack* tr=esd->GetTrack(it);
 243     UInt_t status=tr->GetStatus();
 244     if(!(status&AliESDtrack::kITSrefit)) continue;
 245     if(!(status&AliESDtrack::kTPCin)) continue;
 246     nSelTracks++;
 247   }
 248   fHistoSelTracks->Fill(nSelTracks);
 249
 250   const AliMultiplicity* mult=esd->GetMultiplicity();
 251   Int_t nTracklets=mult->GetNumberOfTracklets();
 252   fHistoTracklets->Fill(nTracklets);
 253
 254   const AliESDVertex *spdv=esd->GetVertex();
 255   if(spdv && spdv->IsFromVertexer3D()){
 256     fHistoSPD3DVtxX->Fill(spdv->GetXv());
 257     fHistoSPD3DVtxY->Fill(spdv->GetYv());
 258     fHistoSPD3DVtxZ->Fill(spdv->GetZv());
 259   }
 260   if(spdv && spdv->IsFromVertexerZ()){
 261     fHistoSPDZVtxX->Fill(spdv->GetXv());
 262     fHistoSPDZVtxY->Fill(spdv->GetYv());
 263     fHistoSPDZVtxZ->Fill(spdv->GetZv());
 264   }
 265   const AliESDVertex *trkv=esd->GetPrimaryVertex();
 266   if(trkv && trkv->GetNContributors()>1){
 267     fHistoTRKVtxX->Fill(trkv->GetXv());
 268     fHistoTRKVtxY->Fill(trkv->GetYv());
 269     fHistoTRKVtxZ->Fill(trkv->GetZv());
 270   }
 271
 272   AliStack* stack=0;
 273
 274   if(fReadMC){
 275     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 276     if (!eventHandler) {
 277       Printf("ERROR: Could not retrieve MC event handler");
 278       return;
 279     }
 280     AliMCEvent* mcEvent = eventHandler->MCEvent();
 281     if (!mcEvent) {
 282       Printf("ERROR: Could not retrieve MC event");
 283       return;
 284     }
 285     stack = mcEvent->Stack();
 286     if (!stack) {
 287       Printf("ERROR: stack not available");
 288       return;
 289     }
 290
 291
 292     Int_t nParticles=stack->GetNtrack();
 293     Double_t dNchdy = 0.;
 294     Int_t nb = 0, nc=0;
 295     Int_t nCharmed=0.;
 296     for (Int_t i=0;i<nParticles;i++){
 297       TParticle* part = (TParticle*)stack->Particle(i);
 298       Int_t absPdg=TMath::Abs(part->GetPdgCode());
 299       if(absPdg==4) nc++;
 300       if(absPdg==5) nb++;
 301       if(stack->IsPhysicalPrimary(i)){
 302         Double_t eta=part->Eta();
 303         if(TMath::Abs(eta)<0.5) dNchdy+=0.6666;   // 2/3 for the ratio charged/all
 304       }
 305       Float_t rapid=-999.;
 306       if (part->Energy() != TMath::Abs(part->Pz())){
 307         rapid=0.5*TMath::Log((part->Energy()+part->Pz())/(part->Energy()-part->Pz()));
 308       }
 309
 310
 311       Int_t iPart=-1;
 312       Int_t iType=0;
 313       Int_t iSpecies=-1;
 314       if(absPdg==421){
 315         iSpecies=0;
 316         iType=CheckD0Decay(i,stack);
 317         if(iType>=0) iPart=0;
 318       }
 319       else if(absPdg==411){
 320         iSpecies=1;
 321         iType=CheckDplusDecay(i,stack);
 322         if(iType>=0) iPart=1;
 323       }
 324       else if(absPdg==413){
 325         iSpecies=2;
 326         iType=CheckDstarDecay(i,stack);
 327         if(iType>=0) iPart=2;
 328       }
 329       else if(absPdg==431){
 330         iSpecies=3;
 331         iType=CheckDsDecay(i,stack);
 332         if(iType==0 || iType==1) iPart=3;
 333       }
 334       else if(absPdg==4122){
 335         iSpecies=4;
 336         iType=CheckLcDecay(i,stack);
 337         if(iType>=0) iPart=4;
 338       }
 339       if(iSpecies<0) continue;
 340       fHistYPtAllDecay[iSpecies]->Fill(part->Pt(),rapid);
 341
 342       TParticle* runningpart=part;
 343       Int_t iFromB=-1;
 344       Int_t pdgmoth=-1;
 345       while(1){
 346         Int_t labmoth=runningpart->GetFirstMother();
 347         if(labmoth==-1) break;
 348         TParticle *mot=(TParticle*)stack->Particle(labmoth);
 349         pdgmoth=TMath::Abs(mot->GetPdgCode());
 350         if(pdgmoth==5){
 351           iFromB=1;
 352           break;
 353         }else if(pdgmoth==4){
 354           iFromB=0;
 355           break;
 356         }
 357         runningpart=mot;
 358       }
 359       if(iFromB<0) continue;
 360       else if(iFromB==0) fHistYPtPromptAllDecay[iSpecies]->Fill(part->Pt(),rapid);
 361       else if(iFromB==1) fHistYPtFeeddownAllDecay[iSpecies]->Fill(part->Pt(),rapid);
 362
 363       if(iPart<0) continue;
 364       if(iType<0) continue;
 365       nCharmed++;
 366       if(iPart==0 && iType<=1){
 367         fHistYPtD0byDecChannel[iType]->Fill(part->Pt(),rapid);
 368       }else if(iPart==1 && iType<=1){
 369         fHistYPtDplusbyDecChannel[iType]->Fill(part->Pt(),rapid);
 370       }else if(iPart==3 && iType<=1){
 371         fHistYPtDsbyDecChannel[iType]->Fill(part->Pt(),rapid);
 372       }
 373
 374       if(iFromB==0 && iPart>=0 && iPart<5) fHistYPtPrompt[iPart]->Fill(part->Pt(),rapid);
 375       else if(iFromB==1 && iPart>=0 && iPart<5) fHistYPtFeeddown[iPart]->Fill(part->Pt(),rapid);
 376
 377     }
 378     fHistoNcharmed->Fill(dNchdy,nCharmed);
 379     fHistoNbVsNc->Fill(nc,nb);
 380   }
 381
 382   PostData(1,fOutput);
 383
 384 }
 385 //______________________________________________________________________________
 386 void AliAnalysisTaskCheckHFMCProd::Terminate(Option_t */*option*/)
 387 {
 388   // Terminate analysis
 389   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 390   if (!fOutput) {
 391     printf("ERROR: fOutput not available\n");
 392     return;
 393   }
 394
 395   return;
 396 }
 397
 398
 399
 400
 401 //______________________________________________________________________________
 402 Int_t AliAnalysisTaskCheckHFMCProd::CheckD0Decay(Int_t labD0, AliStack* stack) const{
 403
 404   if(labD0<0) return -1;
 405   TParticle* dp = (TParticle*)stack->Particle(labD0);
 406   Int_t pdgdp=dp->GetPdgCode();
 407   Int_t nDau=dp->GetNDaughters();
 408
 409   if(nDau==2){
 410     Int_t nKaons=0;
 411     Int_t nPions=0;
 412     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 413       if(iDau<0) return -1;
 414       TParticle* dau=(TParticle*)stack->Particle(iDau);
 415       Int_t pdgdau=dau->GetPdgCode();
 416       if(TMath::Abs(pdgdau)==321){
 417         if(pdgdp>0 && pdgdau>0) return -1;
 418         if(pdgdp<0 && pdgdau<0) return -1;
 419         nKaons++;
 420       }else if(TMath::Abs(pdgdau)==211){
 421         if(pdgdp<0 && pdgdau>0) return -1;
 422         if(pdgdp>0 && pdgdau<0) return -1;
 423         nPions++;
 424       }else{
 425         return -1;
 426       }
 427     }
 428     if(nPions!=1) return -1;
 429     if(nKaons!=1) return -1;
 430     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 431       if(iDau<0) return -1;
 432     }
 433     return 0;
 434   }
 435
 436   if(nDau==3 || nDau==4){
 437     Int_t nKaons=0;
 438     Int_t nPions=0;
 439     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 440       if(iDau<0) return -1;
 441       TParticle* dau=(TParticle*)stack->Particle(iDau);
 442       Int_t pdgdau=dau->GetPdgCode();
 443       if(TMath::Abs(pdgdau)==321){
 444         if(pdgdp>0 && pdgdau>0) return -1;
 445         if(pdgdp<0 && pdgdau<0) return -1;
 446         nKaons++;
 447       }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
 448         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 449           if(resDau<0) return -1;
 450           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 451           Int_t pdgresdau=resdau->GetPdgCode();
 452           if(TMath::Abs(pdgresdau)==321){
 453             if(pdgdp>0 && pdgresdau>0) return -1;
 454             if(pdgdp<0 && pdgresdau<0) return -1;
 455             nKaons++;
 456           }
 457           if(TMath::Abs(pdgresdau)==211){
 458             nPions++;
 459           }
 460         }
 461       }else if(TMath::Abs(pdgdau)==211){
 462           nPions++;
 463       }else{
 464         return -1;
 465       }
 466     }
 467     if(nPions!=3) return -1;
 468     if(nKaons!=1) return -1;
 469     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 470       if(iDau<0) return -1;
 471     }
 472     return 1;
 473   }
 474
 475   return -1;
 476 }
 477
 478
 479 //______________________________________________________________________________
 480 Int_t AliAnalysisTaskCheckHFMCProd::CheckDplusDecay(Int_t labDplus, AliStack* stack) const{
 481
 482   if(labDplus<0) return -1;
 483   TParticle* dp = (TParticle*)stack->Particle(labDplus);
 484   Int_t pdgdp=dp->GetPdgCode();
 485   Int_t nDau=dp->GetNDaughters();
 486
 487   if(nDau==3){
 488     Int_t nKaons=0;
 489     Int_t nPions=0;
 490     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 491       if(iDau<0) return -1;
 492       TParticle* dau=(TParticle*)stack->Particle(iDau);
 493       Int_t pdgdau=dau->GetPdgCode();
 494       if(TMath::Abs(pdgdau)==321){
 495         if(pdgdp>0 && pdgdau>0) return -1;
 496         if(pdgdp<0 && pdgdau<0) return -1;
 497         nKaons++;
 498       }else if(TMath::Abs(pdgdau)==211){
 499         if(pdgdp<0 && pdgdau>0) return -1;
 500         if(pdgdp>0 && pdgdau<0) return -1;
 501         nPions++;
 502       }else{
 503         return -1;
 504       }
 505     }
 506     if(nPions!=2) return -1;
 507     if(nKaons!=1) return -1;
 508     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 509       if(iDau<0) return -1;
 510     }
 511     return 0;
 512   }
 513
 514   if(nDau==2){
 515     Int_t nKaons=0;
 516     Int_t nPions=0;
 517     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 518       if(iDau<0) return -1;
 519       TParticle* dau=(TParticle*)stack->Particle(iDau);
 520       Int_t pdgdau=dau->GetPdgCode();
 521       if(TMath::Abs(pdgdau)==313){
 522         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 523           if(resDau<0) return -1;
 524           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 525           Int_t pdgresdau=resdau->GetPdgCode();
 526           if(TMath::Abs(pdgresdau)==321){
 527             if(pdgdp>0 && pdgresdau>0) return -1;
 528             if(pdgdp<0 && pdgresdau<0) return -1;
 529             nKaons++;
 530           }
 531           if(TMath::Abs(pdgresdau)==211){
 532             if(pdgdp<0 && pdgresdau>0) return -1;
 533             if(pdgdp>0 && pdgresdau<0) return -1;
 534             nPions++;
 535           }
 536         }
 537       }else if(TMath::Abs(pdgdau)==211){
 538           if(pdgdp<0 && pdgdau>0) return -1;
 539           if(pdgdp>0 && pdgdau<0) return -1;
 540           nPions++;
 541       }else{
 542         return -1;
 543       }
 544     }
 545     if(nPions!=2) return -1;
 546     if(nKaons!=1) return -1;
 547     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 548       if(iDau<0) return -1;
 549     }
 550     return 1;
 551   }
 552   return -1;
 553 }
 554
 555 //______________________________________________________________________________
 556 Int_t AliAnalysisTaskCheckHFMCProd::CheckDsDecay(Int_t labDs, AliStack* stack) const{
 557   // Ds decay
 558   if(labDs<0) return -1;
 559   TParticle* dp = (TParticle*)stack->Particle(labDs);
 560   Int_t pdgdp=dp->GetPdgCode();
 561   Int_t nDau=dp->GetNDaughters();
 562
 563   if(nDau==3){
 564     Int_t nKaons=0;
 565     Int_t nPions=0;
 566     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 567       if(iDau<0) return -1;
 568       TParticle* dau=(TParticle*)stack->Particle(iDau);
 569       Int_t pdgdau=dau->GetPdgCode();
 570       if(TMath::Abs(pdgdau)==321){
 571         nKaons++;
 572       }else if(TMath::Abs(pdgdau)==211){
 573         if(pdgdp<0 && pdgdau>0) return -1;
 574         if(pdgdp>0 && pdgdau<0) return -1;
 575         nPions++;
 576       }else{
 577         return -1;
 578       }
 579     }
 580     if(nPions!=1) return -1;
 581     if(nKaons!=2) return -1;
 582     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 583       if(iDau<0) return -1;
 584     }
 585     return 2;
 586   }
 587
 588   if(nDau==2){
 589     Int_t nKaons=0;
 590     Int_t nPions=0;
 591     Bool_t isPhi=kFALSE;
 592     Bool_t isk0st=kFALSE;
 593     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 594       if(iDau<0) return -1;
 595       TParticle* dau=(TParticle*)stack->Particle(iDau);
 596       Int_t pdgdau=dau->GetPdgCode();
 597       if(TMath::Abs(pdgdau)==313){
 598         isk0st=kTRUE;
 599         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 600           if(resDau<0) return -1;
 601           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 602           Int_t pdgresdau=resdau->GetPdgCode();
 603           if(TMath::Abs(pdgresdau)==321){
 604             nKaons++;
 605           }
 606           if(TMath::Abs(pdgresdau)==211){
 607             if(pdgdp<0 && pdgresdau>0) return -1;
 608             if(pdgdp>0 && pdgresdau<0) return -1;
 609             nPions++;
 610           }
 611         }
 612       }else if(TMath::Abs(pdgdau)==333){
 613         isPhi=kTRUE;
 614         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 615           if(resDau<0) return -1;
 616           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 617           if(!resdau) return -1;
 618           Int_t pdgresdau=resdau->GetPdgCode();
 619           if(TMath::Abs(pdgresdau)==321){
 620             nKaons++;
 621           }else{
 622             return -1;
 623           }
 624         }
 625       }else if(TMath::Abs(pdgdau)==211){
 626         if(pdgdp<0 && pdgdau>0) return -1;
 627         if(pdgdp>0 && pdgdau<0) return -1;
 628         nPions++;
 629       }else if(TMath::Abs(pdgdau)==321){
 630         nKaons++;
 631       }else{
 632         return -1;
 633       }
 634     }
 635     if(nPions!=1) return -1;
 636     if(nKaons!=2) return -1;
 637     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 638       if(iDau<0) return -1;
 639     }
 640     if(isk0st) return 1;
 641     else if(isPhi) return 0;
 642     else return 3;
 643   }
 644   return -1;
 645 }
 646
 647 //______________________________________________________________________________
 648 Int_t AliAnalysisTaskCheckHFMCProd::CheckDstarDecay(Int_t labDstar, AliStack* stack) const{
 649
 650   if(labDstar<0) return -1;
 651   TParticle* dp = (TParticle*)stack->Particle(labDstar);
 652   Int_t pdgdp=dp->GetPdgCode();
 653   Int_t nDau=dp->GetNDaughters();
 654   if(nDau!=2) return -1;
 655
 656   Int_t nKaons=0;
 657   Int_t nPions=0;
 658   for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 659     if(iDau<0) return -1;
 660     TParticle* dau=(TParticle*)stack->Particle(iDau);
 661     Int_t pdgdau=dau->GetPdgCode();
 662     if(TMath::Abs(pdgdau)==421){
 663       for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 664         if(resDau<0) return -1;
 665         TParticle* resdau=(TParticle*)stack->Particle(resDau);
 666         Int_t pdgresdau=resdau->GetPdgCode();
 667         if(TMath::Abs(pdgresdau)==321){
 668           if(pdgdp>0 && pdgresdau>0) return -1;
 669           if(pdgdp<0 && pdgresdau<0) return -1;
 670           nKaons++;
 671         }
 672         if(TMath::Abs(pdgresdau)==211){
 673           if(pdgdp<0 && pdgresdau>0) return -1;
 674           if(pdgdp>0 && pdgresdau<0) return -1;
 675           nPions++;
 676         }
 677       }
 678     }else if(TMath::Abs(pdgdau)==211){
 679       if(pdgdp<0 && pdgdau>0) return -1;
 680       if(pdgdp>0 && pdgdau<0) return -1;
 681       nPions++;
 682     }else{
 683       return -1;
 684     }
 685   }
 686   if(nPions!=2) return -1;
 687   if(nKaons!=1) return -1;
 688   for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 689     if(iDau<0) return -1;
 690   }
 691   return 0;
 692
 693 }
 694
 695 //______________________________________________________________________________
 696 Int_t AliAnalysisTaskCheckHFMCProd::CheckLcDecay(Int_t labLc, AliStack* stack) const{
 697   if(labLc<0) return -1;
 698   TParticle* dp = (TParticle*)stack->Particle(labLc);
 699   Int_t pdgdp=dp->GetPdgCode();
 700   Int_t nDau=dp->GetNDaughters();
 701
 702   if(nDau==3){
 703     Int_t nKaons=0;
 704     Int_t nPions=0;
 705     Int_t nProtons=0;
 706     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 707       if(iDau<0) return -1;
 708       TParticle* dau=(TParticle*)stack->Particle(iDau);
 709       Int_t pdgdau=dau->GetPdgCode();
 710       if(TMath::Abs(pdgdau)==321){
 711         if(pdgdp>0 && pdgdau>0) return -1;
 712         if(pdgdp<0 && pdgdau<0) return -1;
 713         nKaons++;
 714       }else if(TMath::Abs(pdgdau)==211){
 715         if(pdgdp<0 && pdgdau>0) return -1;
 716         if(pdgdp>0 && pdgdau<0) return -1;
 717         nPions++;
 718       }else if(TMath::Abs(pdgdau)==2212){
 719         if(pdgdp<0 && pdgdau>0) return -1;
 720         if(pdgdp>0 && pdgdau<0) return -1;
 721         nProtons++;
 722       }else{
 723         return -1;
 724       }
 725     }
 726     if(nPions!=1) return -1;
 727     if(nKaons!=1) return -1;
 728     if(nProtons!=1) return -1;
 729     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 730       if(iDau<0) return -1;
 731     }
 732     return 0;
 733   }
 734
 735   if(nDau==2){
 736     Int_t nKaons=0;
 737     Int_t nPions=0;
 738     Int_t nProtons=0;
 739     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 740       if(iDau<0) return -1;
 741       TParticle* dau=(TParticle*)stack->Particle(iDau);
 742       Int_t pdgdau=dau->GetPdgCode();
 743       if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 || TMath::Abs(pdgdau)==2224
 744          || TMath::Abs(pdgdau)==3122 || TMath::Abs(pdgdau)==311){
 745         Int_t nDauRes=dau->GetNDaughters();
 746         if(nDauRes!=2)  return -1;
 747         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 748           if(resDau<0) return -1;
 749           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 750           Int_t pdgresdau=resdau->GetPdgCode();
 751           if(TMath::Abs(pdgresdau)==321){
 752             if(pdgdp>0 && pdgresdau>0) return -1;
 753             if(pdgdp<0 && pdgresdau<0) return -1;
 754             nKaons++;
 755           }
 756           if(TMath::Abs(pdgresdau)==211){
 757             if(pdgdp<0 && pdgresdau>0) return -1;
 758             if(pdgdp>0 && pdgresdau<0) return -1;
 759             nPions++;
 760           }
 761           if(TMath::Abs(pdgresdau)==2212){
 762             if(pdgdp<0 && pdgresdau>0) return -1;
 763             if(pdgdp>0 && pdgresdau<0) return -1;
 764             nProtons++;
 765           }
 766         }
 767       }else if(TMath::Abs(pdgdau)==321){
 768         if(pdgdp>0 && pdgdau>0) return -1;
 769         if(pdgdp<0 && pdgdau<0) return -1;
 770         nKaons++;
 771       }else if(TMath::Abs(pdgdau)==211){
 772         if(pdgdp<0 && pdgdau>0) return -1;
 773         if(pdgdp>0 && pdgdau<0) return -1;
 774         nPions++;
 775       }else if(TMath::Abs(pdgdau)==2212){
 776         if(pdgdp<0 && pdgdau>0) return -1;
 777         if(pdgdp>0 && pdgdau<0) return -1;
 778         nProtons++;
 779       }else{
 780         return -1;
 781       }
 782     }
 783     if(nPions!=1) return -1;
 784     if(nKaons!=1) return -1;
 785     if(nProtons!=1) return -1;
 786     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 787       if(iDau<0) return -1;
 788     }
 789     return 1;
 790   }
 791   return -1;
 792 }
 793