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   fHistYPtPromptAllD[0] = new TH2F("hyptd0promptAllD","D0 - Prompt",20,0.,20.,40,-2.,2.);
 151   fHistYPtPromptAllD[1] = new TH2F("hyptdpluspromptAllD","Dplus - Prompt",20,0.,20.,40,-2.,2.);
 152   fHistYPtPromptAllD[2] = new TH2F("hyptdstarpromptAllD","Dstar - Prompt",20,0.,20.,40,-2.,2.);
 153   fHistYPtPromptAllD[3] = new TH2F("hyptdspromptAllD","Ds - Prompt",20,0.,20.,40,-2.,2.);
 154   fHistYPtPromptAllD[4] = new TH2F("hyptlcpromptAllD","Lc - Prompt",20,0.,20.,40,-2.,2.);
 155
 156   fHistYPtPrompt[0] = new TH2F("hyptd0prompt","D0 - Prompt",20,0.,20.,20,-2.,2.);
 157   fHistYPtPrompt[1] = new TH2F("hyptdplusprompt","Dplus - Prompt",20,0.,20.,20,-2.,2.);
 158   fHistYPtPrompt[2] = new TH2F("hyptdstarprompt","Dstar - Prompt",20,0.,20.,20,-2.,2.);
 159   fHistYPtPrompt[3] = new TH2F("hyptdsprompt","Ds - Prompt",20,0.,20.,20,-2.,2.);
 160   fHistYPtPrompt[4] = new TH2F("hyptlcprompt","Lc - Prompt",20,0.,20.,20,-2.,2.);
 161
 162   fHistYPtFeeddown[0] = new TH2F("hyptd0feeddown","D0 - Feeddown",20,0.,20.,20,-2.,2.);
 163   fHistYPtFeeddown[1] = new TH2F("hyptdplusfeeddown","Dplus - Feeddown",20,0.,20.,20,-2.,2.);
 164   fHistYPtFeeddown[2] = new TH2F("hyptdstarfeedown","Dstar - Feeddown",20,0.,20.,20,-2.,2.);
 165   fHistYPtFeeddown[3] = new TH2F("hyptdsfeedown","Ds - Feeddown",20,0.,20.,20,-2.,2.);
 166   fHistYPtFeeddown[4] = new TH2F("hyptlcfeedown","Lc - Feeddown",20,0.,20.,20,-2.,2.);
 167
 168   for(Int_t ih=0; ih<5; ih++){
 169     fHistYPtPromptAllD[ih]->Sumw2();
 170     fHistYPtPromptAllD[ih]->SetMinimum(0);
 171     fOutput->Add(fHistYPtPromptAllD[ih]);
 172     fHistYPtPrompt[ih]->Sumw2();
 173     fHistYPtPrompt[ih]->SetMinimum(0);
 174     fOutput->Add(fHistYPtPrompt[ih]);
 175     fHistYPtFeeddown[ih]->Sumw2();
 176     fHistYPtFeeddown[ih]->SetMinimum(0);
 177     fOutput->Add(fHistYPtFeeddown[ih]);
 178   }
 179
 180   fHistYPtD0byDecChannel[0] = new TH2F("hyptD02","D0 - 2prong",20,0.,20.,20,-2.,2.);
 181   fHistYPtD0byDecChannel[1] = new TH2F("hyptD04","D0 - 4prong",20,0.,20.,20,-2.,2.);
 182   fHistYPtDplusbyDecChannel[0] = new TH2F("hyptDplusnonreson","Dplus - non reson",20,0.,20.,20,-2.,2.);
 183   fHistYPtDplusbyDecChannel[1] = new TH2F("hyptDplusreson","Dplus - reson via K0*",20,0.,20.,20,-2.,2.);
 184   fHistYPtDsbyDecChannel[0] = new TH2F("hyptdsphi","Ds - vis Phi",20,0.,20.,20,-2.,2.);
 185   fHistYPtDsbyDecChannel[1] = new TH2F("hyptdsk0st","Ds - via k0*",20,0.,20.,20,-2.,2.);
 186
 187   for(Int_t ih=0; ih<2; ih++){
 188
 189     fHistYPtD0byDecChannel[ih]->Sumw2();
 190     fHistYPtD0byDecChannel[ih]->SetMinimum(0);
 191     fOutput->Add(fHistYPtD0byDecChannel[ih]);
 192     fHistYPtDplusbyDecChannel[ih]->Sumw2();
 193     fHistYPtDplusbyDecChannel[ih]->SetMinimum(0);
 194     fOutput->Add(fHistYPtDplusbyDecChannel[ih]);
 195     fHistYPtDsbyDecChannel[ih]->Sumw2();
 196     fHistYPtDsbyDecChannel[ih]->SetMinimum(0);
 197     fOutput->Add(fHistYPtDsbyDecChannel[ih]);
 198   }
 199
 200
 201   PostData(1,fOutput);
 202
 203 }
 204 //______________________________________________________________________________
 205 void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *)
 206 {
 207   //
 208
 209   AliESDEvent *esd = (AliESDEvent*) (InputEvent());
 210
 211
 212   if(!esd) {
 213     printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
 214     return;
 215   }
 216
 217   fHistoNEvents->Fill(0);
 218
 219   Int_t nTracks=esd->GetNumberOfTracks();
 220   fHistoTracks->Fill(nTracks);
 221   Int_t nSelTracks=0;
 222   for(Int_t it=0; it<nTracks; it++){
 223     AliESDtrack* tr=esd->GetTrack(it);
 224     UInt_t status=tr->GetStatus();
 225     if(!(status&AliESDtrack::kITSrefit)) continue;
 226     if(!(status&AliESDtrack::kTPCin)) continue;
 227     nSelTracks++;
 228   }
 229   fHistoSelTracks->Fill(nSelTracks);
 230
 231   const AliMultiplicity* mult=esd->GetMultiplicity();
 232   Int_t nTracklets=mult->GetNumberOfTracklets();
 233   fHistoTracklets->Fill(nTracklets);
 234
 235   const AliESDVertex *spdv=esd->GetVertex();
 236   if(spdv && spdv->IsFromVertexer3D()){
 237     fHistoSPD3DVtxX->Fill(spdv->GetXv());
 238     fHistoSPD3DVtxY->Fill(spdv->GetYv());
 239     fHistoSPD3DVtxZ->Fill(spdv->GetZv());
 240   }
 241   if(spdv && spdv->IsFromVertexerZ()){
 242     fHistoSPDZVtxX->Fill(spdv->GetXv());
 243     fHistoSPDZVtxY->Fill(spdv->GetYv());
 244     fHistoSPDZVtxZ->Fill(spdv->GetZv());
 245   }
 246   const AliESDVertex *trkv=esd->GetPrimaryVertex();
 247   if(trkv && trkv->GetNContributors()>1){
 248     fHistoTRKVtxX->Fill(trkv->GetXv());
 249     fHistoTRKVtxY->Fill(trkv->GetYv());
 250     fHistoTRKVtxZ->Fill(trkv->GetZv());
 251   }
 252
 253   AliStack* stack=0;
 254
 255   if(fReadMC){
 256     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 257     if (!eventHandler) {
 258       Printf("ERROR: Could not retrieve MC event handler");
 259       return;
 260     }
 261     AliMCEvent* mcEvent = eventHandler->MCEvent();
 262     if (!mcEvent) {
 263       Printf("ERROR: Could not retrieve MC event");
 264       return;
 265     }
 266     stack = mcEvent->Stack();
 267     if (!stack) {
 268       Printf("ERROR: stack not available");
 269       return;
 270     }
 271
 272
 273     Int_t nParticles=stack->GetNtrack();
 274     Double_t dNchdy = 0.;
 275     Int_t nb = 0, nc=0;
 276     Int_t nCharmed=0.;
 277     for (Int_t i=0;i<nParticles;i++){
 278       TParticle* part = (TParticle*)stack->Particle(i);
 279       Int_t absPdg=TMath::Abs(part->GetPdgCode());
 280       if(absPdg==4) nc++;
 281       if(absPdg==5) nb++;
 282       if(stack->IsPhysicalPrimary(i)){
 283         Double_t eta=part->Eta();
 284         if(TMath::Abs(eta)<0.5) dNchdy+=0.6666;   // 2/3 for the ratio charged/all
 285       }
 286       Float_t rapid=-999.;
 287       if (part->Energy() != TMath::Abs(part->Pz())){
 288         rapid=0.5*TMath::Log((part->Energy()+part->Pz())/(part->Energy()-part->Pz()));
 289       }
 290
 291       Int_t iPart=-1;
 292       Int_t iType=0;
 293       if(absPdg==421){
 294         fHistYPtPromptAllD[0]->Fill(part->Pt(),rapid);
 295         iType=CheckD0Decay(i,stack);
 296         if(iType>=0) iPart=0;
 297       }
 298       else if(absPdg==411){
 299         fHistYPtPromptAllD[1]->Fill(part->Pt(),rapid);
 300         iType=CheckDplusDecay(i,stack);
 301         if(iType>=0) iPart=1;
 302       }
 303       else if(absPdg==413){
 304         fHistYPtPromptAllD[2]->Fill(part->Pt(),rapid);
 305         iType=CheckDstarDecay(i,stack);
 306         if(iType>=0) iPart=2;
 307       }
 308       else if(absPdg==431){
 309         fHistYPtPromptAllD[3]->Fill(part->Pt(),rapid);
 310         iType=CheckDsDecay(i,stack);
 311         if(iType==0 || iType==1) iPart=3;
 312       }
 313       else if(absPdg==4122){
 314         fHistYPtPromptAllD[4]->Fill(part->Pt(),rapid);
 315         iType=CheckLcDecay(i,stack);
 316         if(iType>=0) iPart=4;
 317       }
 318       if(iPart<0) continue;
 319       if(iType<0) continue;
 320       nCharmed++;
 321       if(iPart==0 && iType<=1){
 322         fHistYPtD0byDecChannel[iType]->Fill(part->Pt(),rapid);
 323       }else if(iPart==1 && iType<=1){
 324         fHistYPtDplusbyDecChannel[iType]->Fill(part->Pt(),rapid);
 325       }else if(iPart==3 && iType<=1){
 326         fHistYPtDsbyDecChannel[iType]->Fill(part->Pt(),rapid);
 327       }
 328
 329       TParticle* runningpart=part;
 330       Int_t iFromB=-1;
 331       while(1){
 332         Int_t labmoth=runningpart->GetFirstMother();
 333         if(labmoth==-1) break;
 334         TParticle *mot=(TParticle*)stack->Particle(labmoth);
 335         Int_t pdgmoth=TMath::Abs(mot->GetPdgCode());
 336         if(pdgmoth==5){
 337           iFromB=1;
 338           break;
 339         }else if(pdgmoth==4){
 340           iFromB=0;
 341         break;
 342         }
 343         runningpart=mot;
 344       }
 345       if(iFromB<0) continue;
 346       if(iFromB==0 && iPart>=0 && iPart<5) fHistYPtPrompt[iPart]->Fill(part->Pt(),rapid);
 347       else if(iFromB==1 && iPart>=0 && iPart<5) fHistYPtFeeddown[iPart]->Fill(part->Pt(),rapid);
 348
 349     }
 350     printf("  ---> %f %d %d %d\n",dNchdy,nCharmed,nc,nb);
 351     fHistoNcharmed->Fill(dNchdy,nCharmed);
 352     fHistoNbVsNc->Fill(nc,nb);
 353   }
 354
 355   PostData(1,fOutput);
 356
 357 }
 358 //______________________________________________________________________________
 359 void AliAnalysisTaskCheckHFMCProd::Terminate(Option_t */*option*/)
 360 {
 361   // Terminate analysis
 362   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 363   if (!fOutput) {
 364     printf("ERROR: fOutput not available\n");
 365     return;
 366   }
 367
 368   return;
 369 }
 370
 371
 372
 373
 374 //______________________________________________________________________________
 375 Int_t AliAnalysisTaskCheckHFMCProd::CheckD0Decay(Int_t labD0, AliStack* stack) const{
 376
 377   if(labD0<0) return -1;
 378   TParticle* dp = (TParticle*)stack->Particle(labD0);
 379   Int_t pdgdp=dp->GetPdgCode();
 380   Int_t nDau=dp->GetNDaughters();
 381
 382   if(nDau==2){
 383     Int_t nKaons=0;
 384     Int_t nPions=0;
 385     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 386       if(iDau<0) return -1;
 387       TParticle* dau=(TParticle*)stack->Particle(iDau);
 388       Int_t pdgdau=dau->GetPdgCode();
 389       if(TMath::Abs(pdgdau)==321){
 390         if(pdgdp>0 && pdgdau>0) return -1;
 391         if(pdgdp<0 && pdgdau<0) return -1;
 392         nKaons++;
 393       }else if(TMath::Abs(pdgdau)==211){
 394         if(pdgdp<0 && pdgdau>0) return -1;
 395         if(pdgdp>0 && pdgdau<0) return -1;
 396         nPions++;
 397       }else{
 398         return -1;
 399       }
 400     }
 401     if(nPions!=1) return -1;
 402     if(nKaons!=1) return -1;
 403     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 404       if(iDau<0) return -1;
 405     }
 406     return 0;
 407   }
 408
 409   if(nDau==3 || nDau==4){
 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)==113 || TMath::Abs(pdgdau)==313){
 421         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 422           if(resDau<0) return -1;
 423           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 424           Int_t pdgresdau=resdau->GetPdgCode();
 425           if(TMath::Abs(pdgresdau)==321){
 426             if(pdgdp>0 && pdgresdau>0) return -1;
 427             if(pdgdp<0 && pdgresdau<0) return -1;
 428             nKaons++;
 429           }
 430           if(TMath::Abs(pdgresdau)==211){
 431             nPions++;
 432           }
 433         }
 434       }else if(TMath::Abs(pdgdau)==211){
 435           nPions++;
 436       }else{
 437         return -1;
 438       }
 439     }
 440     if(nPions!=3) return -1;
 441     if(nKaons!=1) return -1;
 442     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 443       if(iDau<0) return -1;
 444     }
 445     return 1;
 446   }
 447
 448   return -1;
 449 }
 450
 451
 452 //______________________________________________________________________________
 453 Int_t AliAnalysisTaskCheckHFMCProd::CheckDplusDecay(Int_t labDplus, AliStack* stack) const{
 454
 455   if(labDplus<0) return -1;
 456   TParticle* dp = (TParticle*)stack->Particle(labDplus);
 457   Int_t pdgdp=dp->GetPdgCode();
 458   Int_t nDau=dp->GetNDaughters();
 459
 460   if(nDau==3){
 461     Int_t nKaons=0;
 462     Int_t nPions=0;
 463     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 464       if(iDau<0) return -1;
 465       TParticle* dau=(TParticle*)stack->Particle(iDau);
 466       Int_t pdgdau=dau->GetPdgCode();
 467       if(TMath::Abs(pdgdau)==321){
 468         if(pdgdp>0 && pdgdau>0) return -1;
 469         if(pdgdp<0 && pdgdau<0) return -1;
 470         nKaons++;
 471       }else if(TMath::Abs(pdgdau)==211){
 472         if(pdgdp<0 && pdgdau>0) return -1;
 473         if(pdgdp>0 && pdgdau<0) return -1;
 474         nPions++;
 475       }else{
 476         return -1;
 477       }
 478     }
 479     if(nPions!=2) return -1;
 480     if(nKaons!=1) return -1;
 481     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 482       if(iDau<0) return -1;
 483     }
 484     return 0;
 485   }
 486
 487   if(nDau==2){
 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)==313){
 495         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 496           if(resDau<0) return -1;
 497           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 498           Int_t pdgresdau=resdau->GetPdgCode();
 499           if(TMath::Abs(pdgresdau)==321){
 500             if(pdgdp>0 && pdgresdau>0) return -1;
 501             if(pdgdp<0 && pdgresdau<0) return -1;
 502             nKaons++;
 503           }
 504           if(TMath::Abs(pdgresdau)==211){
 505             if(pdgdp<0 && pdgresdau>0) return -1;
 506             if(pdgdp>0 && pdgresdau<0) return -1;
 507             nPions++;
 508           }
 509         }
 510       }else if(TMath::Abs(pdgdau)==211){
 511           if(pdgdp<0 && pdgdau>0) return -1;
 512           if(pdgdp>0 && pdgdau<0) return -1;
 513           nPions++;
 514       }else{
 515         return -1;
 516       }
 517     }
 518     if(nPions!=2) return -1;
 519     if(nKaons!=1) return -1;
 520     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 521       if(iDau<0) return -1;
 522     }
 523     return 1;
 524   }
 525   return -1;
 526 }
 527
 528 //______________________________________________________________________________
 529 Int_t AliAnalysisTaskCheckHFMCProd::CheckDsDecay(Int_t labDs, AliStack* stack) const{
 530   // Ds decay
 531   if(labDs<0) return -1;
 532   TParticle* dp = (TParticle*)stack->Particle(labDs);
 533   Int_t pdgdp=dp->GetPdgCode();
 534   Int_t nDau=dp->GetNDaughters();
 535
 536   if(nDau==3){
 537     Int_t nKaons=0;
 538     Int_t nPions=0;
 539     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 540       if(iDau<0) return -1;
 541       TParticle* dau=(TParticle*)stack->Particle(iDau);
 542       Int_t pdgdau=dau->GetPdgCode();
 543       if(TMath::Abs(pdgdau)==321){
 544         nKaons++;
 545       }else if(TMath::Abs(pdgdau)==211){
 546         if(pdgdp<0 && pdgdau>0) return -1;
 547         if(pdgdp>0 && pdgdau<0) return -1;
 548         nPions++;
 549       }else{
 550         return -1;
 551       }
 552     }
 553     if(nPions!=1) return -1;
 554     if(nKaons!=2) return -1;
 555     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 556       if(iDau<0) return -1;
 557     }
 558     return 2;
 559   }
 560
 561   if(nDau==2){
 562     Int_t nKaons=0;
 563     Int_t nPions=0;
 564     Bool_t isPhi=kFALSE;
 565     Bool_t isk0st=kFALSE;
 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)==313){
 571         isk0st=kTRUE;
 572         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 573           if(resDau<0) return -1;
 574           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 575           Int_t pdgresdau=resdau->GetPdgCode();
 576           if(TMath::Abs(pdgresdau)==321){
 577             nKaons++;
 578           }
 579           if(TMath::Abs(pdgresdau)==211){
 580             if(pdgdp<0 && pdgresdau>0) return -1;
 581             if(pdgdp>0 && pdgresdau<0) return -1;
 582             nPions++;
 583           }
 584         }
 585       }else if(TMath::Abs(pdgdau)==333){
 586         isPhi=kTRUE;
 587         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 588           if(resDau<0) return -1;
 589           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 590           if(!resdau) return -1;
 591           Int_t pdgresdau=resdau->GetPdgCode();
 592           if(TMath::Abs(pdgresdau)==321){
 593             nKaons++;
 594           }else{
 595             return -1;
 596           }
 597         }
 598       }else if(TMath::Abs(pdgdau)==211){
 599         if(pdgdp<0 && pdgdau>0) return -1;
 600         if(pdgdp>0 && pdgdau<0) return -1;
 601         nPions++;
 602       }else if(TMath::Abs(pdgdau)==321){
 603         nKaons++;
 604       }else{
 605         return -1;
 606       }
 607     }
 608     if(nPions!=1) return -1;
 609     if(nKaons!=2) return -1;
 610     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 611       if(iDau<0) return -1;
 612     }
 613     if(isk0st) return 1;
 614     else if(isPhi) return 0;
 615     else return 3;
 616   }
 617   return -1;
 618 }
 619
 620 //______________________________________________________________________________
 621 Int_t AliAnalysisTaskCheckHFMCProd::CheckDstarDecay(Int_t labDstar, AliStack* stack) const{
 622
 623   if(labDstar<0) return -1;
 624   TParticle* dp = (TParticle*)stack->Particle(labDstar);
 625   Int_t pdgdp=dp->GetPdgCode();
 626   Int_t nDau=dp->GetNDaughters();
 627   if(nDau!=2) return -1;
 628
 629   Int_t nKaons=0;
 630   Int_t nPions=0;
 631   for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 632     if(iDau<0) return -1;
 633     TParticle* dau=(TParticle*)stack->Particle(iDau);
 634     Int_t pdgdau=dau->GetPdgCode();
 635     if(TMath::Abs(pdgdau)==421){
 636       for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 637         if(resDau<0) return -1;
 638         TParticle* resdau=(TParticle*)stack->Particle(resDau);
 639         Int_t pdgresdau=resdau->GetPdgCode();
 640         if(TMath::Abs(pdgresdau)==321){
 641           if(pdgdp>0 && pdgresdau>0) return -1;
 642           if(pdgdp<0 && pdgresdau<0) return -1;
 643           nKaons++;
 644         }
 645         if(TMath::Abs(pdgresdau)==211){
 646           if(pdgdp<0 && pdgresdau>0) return -1;
 647           if(pdgdp>0 && pdgresdau<0) return -1;
 648           nPions++;
 649         }
 650       }
 651     }else if(TMath::Abs(pdgdau)==211){
 652       if(pdgdp<0 && pdgdau>0) return -1;
 653       if(pdgdp>0 && pdgdau<0) return -1;
 654       nPions++;
 655     }else{
 656       return -1;
 657     }
 658   }
 659   if(nPions!=2) return -1;
 660   if(nKaons!=1) return -1;
 661   for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 662     if(iDau<0) return -1;
 663   }
 664   return 0;
 665
 666 }
 667
 668 //______________________________________________________________________________
 669 Int_t AliAnalysisTaskCheckHFMCProd::CheckLcDecay(Int_t labLc, AliStack* stack) const{
 670   if(labLc<0) return -1;
 671   TParticle* dp = (TParticle*)stack->Particle(labLc);
 672   Int_t pdgdp=dp->GetPdgCode();
 673   Int_t nDau=dp->GetNDaughters();
 674
 675   if(nDau==3){
 676     Int_t nKaons=0;
 677     Int_t nPions=0;
 678     Int_t nProtons=0;
 679     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 680       if(iDau<0) return -1;
 681       TParticle* dau=(TParticle*)stack->Particle(iDau);
 682       Int_t pdgdau=dau->GetPdgCode();
 683       if(TMath::Abs(pdgdau)==321){
 684         if(pdgdp>0 && pdgdau>0) return -1;
 685         if(pdgdp<0 && pdgdau<0) return -1;
 686         nKaons++;
 687       }else if(TMath::Abs(pdgdau)==211){
 688         if(pdgdp<0 && pdgdau>0) return -1;
 689         if(pdgdp>0 && pdgdau<0) return -1;
 690         nPions++;
 691       }else if(TMath::Abs(pdgdau)==2212){
 692         if(pdgdp<0 && pdgdau>0) return -1;
 693         if(pdgdp>0 && pdgdau<0) return -1;
 694         nProtons++;
 695       }else{
 696         return -1;
 697       }
 698     }
 699     if(nPions!=1) return -1;
 700     if(nKaons!=1) return -1;
 701     if(nProtons!=1) return -1;
 702     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 703       if(iDau<0) return -1;
 704     }
 705     return 0;
 706   }
 707
 708   if(nDau==2){
 709     Int_t nKaons=0;
 710     Int_t nPions=0;
 711     Int_t nProtons=0;
 712     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 713       if(iDau<0) return -1;
 714       TParticle* dau=(TParticle*)stack->Particle(iDau);
 715       Int_t pdgdau=dau->GetPdgCode();
 716       if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 || TMath::Abs(pdgdau)==2224
 717          || TMath::Abs(pdgdau)==3122 || TMath::Abs(pdgdau)==311){
 718         Int_t nDauRes=dau->GetNDaughters();
 719         if(nDauRes!=2)  return -1;
 720         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 721           if(resDau<0) return -1;
 722           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 723           Int_t pdgresdau=resdau->GetPdgCode();
 724           if(TMath::Abs(pdgresdau)==321){
 725             if(pdgdp>0 && pdgresdau>0) return -1;
 726             if(pdgdp<0 && pdgresdau<0) return -1;
 727             nKaons++;
 728           }
 729           if(TMath::Abs(pdgresdau)==211){
 730             if(pdgdp<0 && pdgresdau>0) return -1;
 731             if(pdgdp>0 && pdgresdau<0) return -1;
 732             nPions++;
 733           }
 734           if(TMath::Abs(pdgresdau)==2212){
 735             if(pdgdp<0 && pdgresdau>0) return -1;
 736             if(pdgdp>0 && pdgresdau<0) return -1;
 737             nProtons++;
 738           }
 739         }
 740       }else if(TMath::Abs(pdgdau)==321){
 741         if(pdgdp>0 && pdgdau>0) return -1;
 742         if(pdgdp<0 && pdgdau<0) return -1;
 743         nKaons++;
 744       }else if(TMath::Abs(pdgdau)==211){
 745         if(pdgdp<0 && pdgdau>0) return -1;
 746         if(pdgdp>0 && pdgdau<0) return -1;
 747         nPions++;
 748       }else if(TMath::Abs(pdgdau)==2212){
 749         if(pdgdp<0 && pdgdau>0) return -1;
 750         if(pdgdp>0 && pdgdau<0) return -1;
 751         nProtons++;
 752       }else{
 753         return -1;
 754       }
 755     }
 756     if(nPions!=1) return -1;
 757     if(nKaons!=1) return -1;
 758     if(nProtons!=1) return -1;
 759     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 760       if(iDau<0) return -1;
 761     }
 762     return 1;
 763   }
 764   return -1;
 765 }
 766