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   fHistoPhysPrim(0),
  52   fHistoTracks(0),
  53   fHistoSelTracks(0),
  54   fHistoTracklets(0),
  55   fHistoSPD3DVtxX(0),
  56   fHistoSPD3DVtxY(0),
  57   fHistoSPD3DVtxZ(0),
  58   fHistoSPDZVtxX(0),
  59   fHistoSPDZVtxY(0),
  60   fHistoSPDZVtxZ(0),
  61   fHistoTRKVtxX(0),
  62   fHistoTRKVtxY(0),
  63   fHistoTRKVtxZ(0),
  64   fHistoNcharmed(0),
  65   fHistoNbVsNc(0),
  66   fHistOriginPrompt(0),
  67   fHistOriginFeeddown(0),
  68   fHistMotherID(0),
  69   fSearchUpToQuark(kFALSE),
  70   fSystem(0),
  71   fReadMC(kTRUE)
  72 {
  73   //
  74   DefineInput(0, TChain::Class());
  75   DefineOutput(1, TList::Class());
  76 }
  77
  78
  79 //___________________________________________________________________________
  80 AliAnalysisTaskCheckHFMCProd::~AliAnalysisTaskCheckHFMCProd(){
  81   //
  82   if (AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return;
  83   if (fOutput) {
  84     delete fOutput;
  85     fOutput = 0;
  86   }
  87 }
  88
  89 //___________________________________________________________________________
  90 void AliAnalysisTaskCheckHFMCProd::UserCreateOutputObjects() {
  91   // create output histos
  92
  93   fOutput = new TList();
  94   fOutput->SetOwner();
  95   fOutput->SetName("OutputHistos");
  96
  97   fHistoNEvents = new TH1F("hNEvents", "Number of processed events",3,-0.5,2.5);
  98   fHistoNEvents->Sumw2();
  99   fHistoNEvents->SetMinimum(0);
 100   fOutput->Add(fHistoNEvents);
 101
 102   Double_t maxMult=100.;
 103   if(fSystem==1) maxMult=10000.;
 104   if(fSystem==2) maxMult=500.;
 105   fHistoPhysPrim = new TH1F("hPhysPrim","",100,0.,maxMult);
 106   fHistoPhysPrim->Sumw2();
 107   fOutput->Add(fHistoPhysPrim);
 108   fHistoTracks = new TH1F("hTracks","",100,0.,maxMult*2);
 109   fHistoTracks->Sumw2();
 110   fOutput->Add(fHistoTracks);
 111   fHistoSelTracks = new TH1F("hSelTracks","",100,0.,maxMult);
 112   fHistoSelTracks->Sumw2();
 113   fOutput->Add(fHistoSelTracks);
 114   fHistoTracklets = new TH1F("hTracklets","",100,0.,maxMult);
 115   fHistoTracklets->Sumw2();
 116   fOutput->Add(fHistoTracklets);
 117
 118   fHistoSPD3DVtxX = new TH1F("hSPD3DvX","",100,-1.,1.);
 119   fHistoSPD3DVtxX->Sumw2();
 120   fOutput->Add(fHistoSPD3DVtxX);
 121   fHistoSPD3DVtxY = new TH1F("hSPD3DvY","",100,-1.,1.);
 122   fHistoSPD3DVtxY->Sumw2();
 123   fOutput->Add(fHistoSPD3DVtxY);
 124   fHistoSPD3DVtxZ = new TH1F("hSPD3DvZ","",100,-15.,15.);
 125   fHistoSPD3DVtxZ->Sumw2();
 126   fOutput->Add(fHistoSPD3DVtxZ);
 127
 128   fHistoSPDZVtxX = new TH1F("hSPDZvX","",100,-1.,1.);
 129   fHistoSPDZVtxX->Sumw2();
 130   fOutput->Add(fHistoSPDZVtxX);
 131   fHistoSPDZVtxY = new TH1F("hSPDZvY","",100,-1.,1.);
 132   fHistoSPDZVtxY->Sumw2();
 133   fOutput->Add(fHistoSPDZVtxY);
 134   fHistoSPDZVtxZ = new TH1F("hSPDZvZ","",100,-15.,15.);
 135   fHistoSPDZVtxZ->Sumw2();
 136   fOutput->Add(fHistoSPDZVtxZ);
 137
 138
 139   fHistoTRKVtxX = new TH1F("hTRKvX","",100,-1.,1.);
 140   fHistoTRKVtxX->Sumw2();
 141   fOutput->Add(fHistoTRKVtxX);
 142   fHistoTRKVtxY = new TH1F("hTRKvY","",100,-1.,1.);
 143   fHistoTRKVtxY->Sumw2();
 144   fOutput->Add(fHistoTRKVtxY);
 145   fHistoTRKVtxZ = new TH1F("hTRKvZ","",100,-15.,15.);
 146   fHistoTRKVtxZ->Sumw2();
 147   fOutput->Add(fHistoTRKVtxZ);
 148
 149   Int_t nBinscb=11;
 150   if(fSystem==1) nBinscb=400;
 151   if(fSystem==2) nBinscb=21;
 152   Double_t maxncn=nBinscb-0.5;
 153   fHistoNcharmed = new TH2F("hncharmed","",100,0.,maxMult,nBinscb,-0.5,maxncn);
 154   fHistoNcharmed->Sumw2();
 155   fOutput->Add(fHistoNcharmed);
 156   fHistoNbVsNc = new TH2F("hnbvsnc","",nBinscb,-0.5,maxncn,nBinscb,-0.5,maxncn);
 157   fHistoNbVsNc->Sumw2();
 158   fOutput->Add(fHistoNbVsNc);
 159
 160   fHistYPtPrompt[0] = new TH2F("hyptD0prompt","D0 - Prompt",40,0.,40.,20,-2.,2.);
 161   fHistYPtPrompt[1] = new TH2F("hyptDplusprompt","Dplus - Prompt",40,0.,40.,20,-2.,2.);
 162   fHistYPtPrompt[2] = new TH2F("hyptDstarprompt","Dstar - Prompt",40,0.,40.,20,-2.,2.);
 163   fHistYPtPrompt[3] = new TH2F("hyptDsprompt","Ds - Prompt",40,0.,40.,20,-2.,2.);
 164   fHistYPtPrompt[4] = new TH2F("hyptLcprompt","Lc - Prompt",40,0.,40.,20,-2.,2.);
 165
 166   fHistBYPtAllDecay[0] = new TH2F("hyptB0AllDecay","B0 - All",40,0.,40.,40,-2.,2.);
 167   fHistBYPtAllDecay[1] = new TH2F("hyptBplusAllDecay","Bplus - All",40,0.,40.,40,-2.,2.);
 168   fHistBYPtAllDecay[2] = new TH2F("hyptBstarAllDecay","Bstar - All",40,0.,40.,40,-2.,2.);
 169   fHistBYPtAllDecay[3] = new TH2F("hyptBsAllDecay","Bs - All",40,0.,40.,40,-2.,2.);
 170   fHistBYPtAllDecay[4] = new TH2F("hyptLbAllDecay","LB - All",40,0.,40.,40,-2.,2.);
 171
 172   fHistYPtAllDecay[0] = new TH2F("hyptD0AllDecay","D0 - All",40,0.,40.,40,-2.,2.);
 173   fHistYPtAllDecay[1] = new TH2F("hyptDplusAllDecay","Dplus - All",40,0.,40.,40,-2.,2.);
 174   fHistYPtAllDecay[2] = new TH2F("hyptDstarAllDecay","Dstar - All",40,0.,40.,40,-2.,2.);
 175   fHistYPtAllDecay[3] = new TH2F("hyptDsAllDecay","Ds - All",40,0.,40.,40,-2.,2.);
 176   fHistYPtAllDecay[4] = new TH2F("hyptLcAllDecay","Lc - All",40,0.,40.,40,-2.,2.);
 177
 178   fHistYPtPromptAllDecay[0] = new TH2F("hyptD0promptAllDecay","D0 - Prompt",40,0.,40.,40,-2.,2.);
 179   fHistYPtPromptAllDecay[1] = new TH2F("hyptDpluspromptAllDecay","Dplus - Prompt",40,0.,40.,40,-2.,2.);
 180   fHistYPtPromptAllDecay[2] = new TH2F("hyptDstarpromptAllDecay","Dstar - Prompt",40,0.,40.,40,-2.,2.);
 181   fHistYPtPromptAllDecay[3] = new TH2F("hyptDspromptAllDecay","Ds - Prompt",40,0.,40.,40,-2.,2.);
 182   fHistYPtPromptAllDecay[4] = new TH2F("hyptLcpromptAllDecay","Lc - Prompt",40,0.,40.,40,-2.,2.);
 183
 184   fHistYPtFeeddownAllDecay[0] = new TH2F("hyptD0feeddownAllDecay","D0 - FromB",40,0.,40.,40,-2.,2.);
 185   fHistYPtFeeddownAllDecay[1] = new TH2F("hyptDplusfeeddownAllDecay","Dplus - FromB",40,0.,40.,40,-2.,2.);
 186   fHistYPtFeeddownAllDecay[2] = new TH2F("hyptDstarfeeddownAllDecay","Dstar - FromB",40,0.,40.,40,-2.,2.);
 187   fHistYPtFeeddownAllDecay[3] = new TH2F("hyptDsfeeddownAllDecay","Ds - FromB",40,0.,40.,40,-2.,2.);
 188   fHistYPtFeeddownAllDecay[4] = new TH2F("hyptLcfeeddownAllDecay","Lc - FromB",40,0.,40.,40,-2.,2.);
 189
 190
 191   fHistYPtFeeddown[0] = new TH2F("hyptD0feeddown","D0 - Feeddown",40,0.,40.,20,-2.,2.);
 192   fHistYPtFeeddown[1] = new TH2F("hyptDplusfeeddown","Dplus - Feeddown",40,0.,40.,20,-2.,2.);
 193   fHistYPtFeeddown[2] = new TH2F("hyptDstarfeedown","Dstar - Feeddown",40,0.,40.,20,-2.,2.);
 194   fHistYPtFeeddown[3] = new TH2F("hyptDsfeedown","Ds - Feeddown",40,0.,40.,20,-2.,2.);
 195   fHistYPtFeeddown[4] = new TH2F("hyptLcfeedown","Lc - Feeddown",40,0.,40.,20,-2.,2.);
 196
 197   for(Int_t ih=0; ih<5; ih++){
 198     fHistBYPtAllDecay[ih]->Sumw2();
 199     fHistBYPtAllDecay[ih]->SetMinimum(0);
 200     fOutput->Add(fHistBYPtAllDecay[ih]);
 201     fHistYPtAllDecay[ih]->Sumw2();
 202     fHistYPtAllDecay[ih]->SetMinimum(0);
 203     fOutput->Add(fHistYPtAllDecay[ih]);
 204     fHistYPtPromptAllDecay[ih]->Sumw2();
 205     fHistYPtPromptAllDecay[ih]->SetMinimum(0);
 206     fOutput->Add(fHistYPtPromptAllDecay[ih]);
 207     fHistYPtFeeddownAllDecay[ih]->Sumw2();
 208     fHistYPtFeeddownAllDecay[ih]->SetMinimum(0);
 209     fOutput->Add(fHistYPtFeeddownAllDecay[ih]);
 210     fHistYPtPrompt[ih]->Sumw2();
 211     fHistYPtPrompt[ih]->SetMinimum(0);
 212     fOutput->Add(fHistYPtPrompt[ih]);
 213     fHistYPtFeeddown[ih]->Sumw2();
 214     fHistYPtFeeddown[ih]->SetMinimum(0);
 215     fOutput->Add(fHistYPtFeeddown[ih]);
 216   }
 217
 218   fHistYPtD0byDecChannel[0] = new TH2F("hyptD02","D0 - 2prong",40,0.,40.,20,-2.,2.);
 219   fHistYPtD0byDecChannel[1] = new TH2F("hyptD04","D0 - 4prong",40,0.,40.,20,-2.,2.);
 220   fHistYPtDplusbyDecChannel[0] = new TH2F("hyptDplusnonreson","Dplus - non reson",40,0.,40.,20,-2.,2.);
 221   fHistYPtDplusbyDecChannel[1] = new TH2F("hyptDplusreson","Dplus - reson via K0*",40,0.,40.,20,-2.,2.);
 222   fHistYPtDsbyDecChannel[0] = new TH2F("hyptDsphi","Ds - vis Phi",40,0.,40.,20,-2.,2.);
 223   fHistYPtDsbyDecChannel[1] = new TH2F("hyptDsk0st","Ds - via k0*",40,0.,40.,20,-2.,2.);
 224
 225   for(Int_t ih=0; ih<2; ih++){
 226
 227     fHistYPtD0byDecChannel[ih]->Sumw2();
 228     fHistYPtD0byDecChannel[ih]->SetMinimum(0);
 229     fOutput->Add(fHistYPtD0byDecChannel[ih]);
 230     fHistYPtDplusbyDecChannel[ih]->Sumw2();
 231     fHistYPtDplusbyDecChannel[ih]->SetMinimum(0);
 232     fOutput->Add(fHistYPtDplusbyDecChannel[ih]);
 233     fHistYPtDsbyDecChannel[ih]->Sumw2();
 234     fHistYPtDsbyDecChannel[ih]->SetMinimum(0);
 235     fOutput->Add(fHistYPtDsbyDecChannel[ih]);
 236   }
 237
 238   fHistOriginPrompt=new TH1F("hOriginPrompt","",100,0.,0.5);
 239   fHistOriginPrompt->Sumw2();
 240   fHistOriginPrompt->SetMinimum(0);
 241   fOutput->Add(fHistOriginPrompt);
 242   fHistOriginFeeddown=new TH1F("hOriginFeeddown","",100,0.,0.5);
 243   fHistOriginFeeddown->Sumw2();
 244   fHistOriginFeeddown->SetMinimum(0);
 245   fOutput->Add(fHistOriginFeeddown);
 246   fHistMotherID=new TH1F("hMotherID","",1000,-1.5,998.5);
 247   fHistMotherID->SetMinimum(0);
 248   fOutput->Add(fHistMotherID);
 249
 250   PostData(1,fOutput);
 251
 252 }
 253 //______________________________________________________________________________
 254 void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *)
 255 {
 256   //
 257
 258   AliESDEvent *esd = (AliESDEvent*) (InputEvent());
 259
 260
 261   if(!esd) {
 262     printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
 263     return;
 264   }
 265
 266   fHistoNEvents->Fill(0);
 267
 268   Int_t nTracks=esd->GetNumberOfTracks();
 269   fHistoTracks->Fill(nTracks);
 270   Int_t nSelTracks=0;
 271   for(Int_t it=0; it<nTracks; it++){
 272     AliESDtrack* tr=esd->GetTrack(it);
 273     UInt_t status=tr->GetStatus();
 274     if(!(status&AliESDtrack::kITSrefit)) continue;
 275     if(!(status&AliESDtrack::kTPCin)) continue;
 276     nSelTracks++;
 277   }
 278   fHistoSelTracks->Fill(nSelTracks);
 279
 280   const AliMultiplicity* mult=esd->GetMultiplicity();
 281   Int_t nTracklets=mult->GetNumberOfTracklets();
 282   fHistoTracklets->Fill(nTracklets);
 283
 284   const AliESDVertex *spdv=esd->GetVertex();
 285   if(spdv && spdv->IsFromVertexer3D()){
 286     fHistoSPD3DVtxX->Fill(spdv->GetXv());
 287     fHistoSPD3DVtxY->Fill(spdv->GetYv());
 288     fHistoSPD3DVtxZ->Fill(spdv->GetZv());
 289   }
 290   if(spdv && spdv->IsFromVertexerZ()){
 291     fHistoSPDZVtxX->Fill(spdv->GetXv());
 292     fHistoSPDZVtxY->Fill(spdv->GetYv());
 293     fHistoSPDZVtxZ->Fill(spdv->GetZv());
 294   }
 295   const AliESDVertex *trkv=esd->GetPrimaryVertex();
 296   if(trkv && trkv->GetNContributors()>1){
 297     fHistoTRKVtxX->Fill(trkv->GetXv());
 298     fHistoTRKVtxY->Fill(trkv->GetYv());
 299     fHistoTRKVtxZ->Fill(trkv->GetZv());
 300   }
 301
 302   AliStack* stack=0;
 303
 304   if(fReadMC){
 305     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 306     if (!eventHandler) {
 307       Printf("ERROR: Could not retrieve MC event handler");
 308       return;
 309     }
 310     AliMCEvent* mcEvent = eventHandler->MCEvent();
 311     if (!mcEvent) {
 312       Printf("ERROR: Could not retrieve MC event");
 313       return;
 314     }
 315     stack = mcEvent->Stack();
 316     if (!stack) {
 317       Printf("ERROR: stack not available");
 318       return;
 319     }
 320     const AliVVertex* mcVert=mcEvent->GetPrimaryVertex();
 321     if(!mcVert){
 322       Printf("ERROR: generated vertex not available");
 323       return;
 324     }
 325
 326
 327     Int_t nParticles=stack->GetNtrack();
 328     Double_t dNchdy = 0.;
 329     Int_t nb = 0, nc=0;
 330     Int_t nCharmed=0;
 331     Int_t nPhysPrim=0;
 332     for (Int_t i=0;i<nParticles;i++){
 333       TParticle* part = (TParticle*)stack->Particle(i);
 334       Int_t absPdg=TMath::Abs(part->GetPdgCode());
 335       if(absPdg==4) nc++;
 336       if(absPdg==5) nb++;
 337       if(stack->IsPhysicalPrimary(i)){
 338         Double_t eta=part->Eta();
 339         if(TMath::Abs(eta)<0.5){
 340           dNchdy+=0.6666;   // 2/3 for the ratio charged/all
 341           nPhysPrim++;
 342         }
 343       }
 344       Float_t rapid=-999.;
 345       if (part->Energy() != TMath::Abs(part->Pz())){
 346         rapid=0.5*TMath::Log((part->Energy()+part->Pz())/(part->Energy()-part->Pz()));
 347       }
 348
 349
 350       Int_t iPart=-1;
 351       Int_t iType=0;
 352       Int_t iSpecies=-1;
 353       if(absPdg==421){
 354         iSpecies=0;
 355         iType=CheckD0Decay(i,stack);
 356         if(iType>=0) iPart=0;
 357       }
 358       else if(absPdg==411){
 359         iSpecies=1;
 360         iType=CheckDplusDecay(i,stack);
 361         if(iType>=0) iPart=1;
 362       }
 363       else if(absPdg==413){
 364         iSpecies=2;
 365         iType=CheckDstarDecay(i,stack);
 366         if(iType>=0) iPart=2;
 367       }
 368       else if(absPdg==431){
 369         iSpecies=3;
 370         iType=CheckDsDecay(i,stack);
 371         if(iType==0 || iType==1) iPart=3;
 372       }
 373       else if(absPdg==4122){
 374         iSpecies=4;
 375         iType=CheckLcDecay(i,stack);
 376         if(iType>=0) iPart=4;
 377       }
 378       if(iSpecies>=0) fHistYPtAllDecay[iSpecies]->Fill(part->Pt(),rapid);
 379
 380       // check beauty mesons
 381       if(absPdg==511) fHistBYPtAllDecay[0]->Fill(part->Pt(),rapid);
 382       else if(absPdg==521) fHistBYPtAllDecay[1]->Fill(part->Pt(),rapid);
 383       else if(absPdg==513) fHistBYPtAllDecay[2]->Fill(part->Pt(),rapid);
 384       else if(absPdg==531) fHistBYPtAllDecay[3]->Fill(part->Pt(),rapid);
 385       else if(absPdg==5122) fHistBYPtAllDecay[4]->Fill(part->Pt(),rapid);
 386
 387       if(iSpecies<0) continue; // not a charm meson
 388
 389       Double_t distx=part->Vx()-mcVert->GetX();
 390       Double_t disty=part->Vy()-mcVert->GetY();
 391       Double_t distz=part->Vz()-mcVert->GetZ();
 392       Double_t distToVert=TMath::Sqrt(distx*distx+disty*disty+distz*distz);
 393       fHistMotherID->Fill(part->GetFirstMother());
 394       TParticle* runningpart=part;
 395       Int_t iFromB=-1;
 396       Int_t pdgmoth=-1;
 397       if(fSearchUpToQuark){
 398         while(1){
 399           Int_t labmoth=runningpart->GetFirstMother();
 400           if(labmoth==-1) break;
 401           TParticle *mot=(TParticle*)stack->Particle(labmoth);
 402           pdgmoth=TMath::Abs(mot->GetPdgCode());
 403           if(pdgmoth==5){
 404             iFromB=1;
 405             break;
 406           }else if(pdgmoth==4){
 407             iFromB=0;
 408             break;
 409           }
 410           runningpart=mot;
 411         }
 412       }else{
 413         iFromB=0;
 414         while(1){
 415           Int_t labmoth=runningpart->GetFirstMother();
 416           if(labmoth==-1) break;
 417           TParticle *mot=(TParticle*)stack->Particle(labmoth);
 418           pdgmoth=TMath::Abs(mot->GetPdgCode());
 419           if(pdgmoth>=500 && pdgmoth<=599){
 420             iFromB=1;
 421             break;
 422           }
 423           if(pdgmoth>=5000 && pdgmoth<=5999){
 424             iFromB=1;
 425             break;
 426           }
 427           runningpart=mot;
 428         }
 429       }
 430
 431       if(iFromB==0){
 432         fHistYPtPromptAllDecay[iSpecies]->Fill(part->Pt(),rapid);
 433         fHistOriginPrompt->Fill(distToVert);
 434       }
 435       else if(iFromB==1){
 436         fHistYPtFeeddownAllDecay[iSpecies]->Fill(part->Pt(),rapid);
 437         fHistOriginFeeddown->Fill(distToVert);
 438       }
 439
 440       if(iPart<0) continue;
 441       if(iType<0) continue;
 442       nCharmed++;
 443       if(iPart==0 && iType<=1){
 444         fHistYPtD0byDecChannel[iType]->Fill(part->Pt(),rapid);
 445       }else if(iPart==1 && iType<=1){
 446         fHistYPtDplusbyDecChannel[iType]->Fill(part->Pt(),rapid);
 447       }else if(iPart==3 && iType<=1){
 448         fHistYPtDsbyDecChannel[iType]->Fill(part->Pt(),rapid);
 449       }
 450
 451       if(iFromB==0 && iPart>=0 && iPart<5) fHistYPtPrompt[iPart]->Fill(part->Pt(),rapid);
 452       else if(iFromB==1 && iPart>=0 && iPart<5) fHistYPtFeeddown[iPart]->Fill(part->Pt(),rapid);
 453     }
 454
 455     fHistoNcharmed->Fill(dNchdy,nCharmed);
 456     fHistoNbVsNc->Fill(nc,nb);
 457     fHistoPhysPrim->Fill(nPhysPrim);
 458   }
 459
 460   PostData(1,fOutput);
 461
 462 }
 463 //______________________________________________________________________________
 464 void AliAnalysisTaskCheckHFMCProd::Terminate(Option_t */*option*/)
 465 {
 466   // Terminate analysis
 467   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 468   if (!fOutput) {
 469     printf("ERROR: fOutput not available\n");
 470     return;
 471   }
 472
 473   return;
 474 }
 475
 476
 477
 478
 479 //______________________________________________________________________________
 480 Int_t AliAnalysisTaskCheckHFMCProd::CheckD0Decay(Int_t labD0, AliStack* stack) const{
 481
 482   if(labD0<0) return -1;
 483   TParticle* dp = (TParticle*)stack->Particle(labD0);
 484   Int_t pdgdp=dp->GetPdgCode();
 485   Int_t nDau=dp->GetNDaughters();
 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)==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!=1) 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==3 || nDau==4){
 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)==321){
 522         if(pdgdp>0 && pdgdau>0) return -1;
 523         if(pdgdp<0 && pdgdau<0) return -1;
 524         nKaons++;
 525       }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
 526         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 527           if(resDau<0) return -1;
 528           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 529           Int_t pdgresdau=resdau->GetPdgCode();
 530           if(TMath::Abs(pdgresdau)==321){
 531             if(pdgdp>0 && pdgresdau>0) return -1;
 532             if(pdgdp<0 && pdgresdau<0) return -1;
 533             nKaons++;
 534           }
 535           if(TMath::Abs(pdgresdau)==211){
 536             nPions++;
 537           }
 538         }
 539       }else if(TMath::Abs(pdgdau)==211){
 540           nPions++;
 541       }else{
 542         return -1;
 543       }
 544     }
 545     if(nPions!=3) 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
 553   return -1;
 554 }
 555
 556
 557 //______________________________________________________________________________
 558 Int_t AliAnalysisTaskCheckHFMCProd::CheckDplusDecay(Int_t labDplus, AliStack* stack) const{
 559
 560   if(labDplus<0) return -1;
 561   TParticle* dp = (TParticle*)stack->Particle(labDplus);
 562   Int_t pdgdp=dp->GetPdgCode();
 563   Int_t nDau=dp->GetNDaughters();
 564
 565   if(nDau==3){
 566     Int_t nKaons=0;
 567     Int_t nPions=0;
 568     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 569       if(iDau<0) return -1;
 570       TParticle* dau=(TParticle*)stack->Particle(iDau);
 571       Int_t pdgdau=dau->GetPdgCode();
 572       if(TMath::Abs(pdgdau)==321){
 573         if(pdgdp>0 && pdgdau>0) return -1;
 574         if(pdgdp<0 && pdgdau<0) return -1;
 575         nKaons++;
 576       }else if(TMath::Abs(pdgdau)==211){
 577         if(pdgdp<0 && pdgdau>0) return -1;
 578         if(pdgdp>0 && pdgdau<0) return -1;
 579         nPions++;
 580       }else{
 581         return -1;
 582       }
 583     }
 584     if(nPions!=2) return -1;
 585     if(nKaons!=1) return -1;
 586     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 587       if(iDau<0) return -1;
 588     }
 589     return 0;
 590   }
 591
 592   if(nDau==2){
 593     Int_t nKaons=0;
 594     Int_t nPions=0;
 595     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 596       if(iDau<0) return -1;
 597       TParticle* dau=(TParticle*)stack->Particle(iDau);
 598       Int_t pdgdau=dau->GetPdgCode();
 599       if(TMath::Abs(pdgdau)==313){
 600         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 601           if(resDau<0) return -1;
 602           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 603           Int_t pdgresdau=resdau->GetPdgCode();
 604           if(TMath::Abs(pdgresdau)==321){
 605             if(pdgdp>0 && pdgresdau>0) return -1;
 606             if(pdgdp<0 && pdgresdau<0) return -1;
 607             nKaons++;
 608           }
 609           if(TMath::Abs(pdgresdau)==211){
 610             if(pdgdp<0 && pdgresdau>0) return -1;
 611             if(pdgdp>0 && pdgresdau<0) return -1;
 612             nPions++;
 613           }
 614         }
 615       }else if(TMath::Abs(pdgdau)==211){
 616           if(pdgdp<0 && pdgdau>0) return -1;
 617           if(pdgdp>0 && pdgdau<0) return -1;
 618           nPions++;
 619       }else{
 620         return -1;
 621       }
 622     }
 623     if(nPions!=2) return -1;
 624     if(nKaons!=1) return -1;
 625     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 626       if(iDau<0) return -1;
 627     }
 628     return 1;
 629   }
 630   return -1;
 631 }
 632
 633 //______________________________________________________________________________
 634 Int_t AliAnalysisTaskCheckHFMCProd::CheckDsDecay(Int_t labDs, AliStack* stack) const{
 635   // Ds decay
 636   if(labDs<0) return -1;
 637   TParticle* dp = (TParticle*)stack->Particle(labDs);
 638   Int_t pdgdp=dp->GetPdgCode();
 639   Int_t nDau=dp->GetNDaughters();
 640
 641   if(nDau==3){
 642     Int_t nKaons=0;
 643     Int_t nPions=0;
 644     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 645       if(iDau<0) return -1;
 646       TParticle* dau=(TParticle*)stack->Particle(iDau);
 647       Int_t pdgdau=dau->GetPdgCode();
 648       if(TMath::Abs(pdgdau)==321){
 649         nKaons++;
 650       }else if(TMath::Abs(pdgdau)==211){
 651         if(pdgdp<0 && pdgdau>0) return -1;
 652         if(pdgdp>0 && pdgdau<0) return -1;
 653         nPions++;
 654       }else{
 655         return -1;
 656       }
 657     }
 658     if(nPions!=1) return -1;
 659     if(nKaons!=2) return -1;
 660     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 661       if(iDau<0) return -1;
 662     }
 663     return 2;
 664   }
 665
 666   if(nDau==2){
 667     Int_t nKaons=0;
 668     Int_t nPions=0;
 669     Bool_t isPhi=kFALSE;
 670     Bool_t isk0st=kFALSE;
 671     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 672       if(iDau<0) return -1;
 673       TParticle* dau=(TParticle*)stack->Particle(iDau);
 674       Int_t pdgdau=dau->GetPdgCode();
 675       if(TMath::Abs(pdgdau)==313){
 676         isk0st=kTRUE;
 677         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 678           if(resDau<0) return -1;
 679           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 680           Int_t pdgresdau=resdau->GetPdgCode();
 681           if(TMath::Abs(pdgresdau)==321){
 682             nKaons++;
 683           }
 684           if(TMath::Abs(pdgresdau)==211){
 685             if(pdgdp<0 && pdgresdau>0) return -1;
 686             if(pdgdp>0 && pdgresdau<0) return -1;
 687             nPions++;
 688           }
 689         }
 690       }else if(TMath::Abs(pdgdau)==333){
 691         isPhi=kTRUE;
 692         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 693           if(resDau<0) return -1;
 694           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 695           if(!resdau) return -1;
 696           Int_t pdgresdau=resdau->GetPdgCode();
 697           if(TMath::Abs(pdgresdau)==321){
 698             nKaons++;
 699           }else{
 700             return -1;
 701           }
 702         }
 703       }else if(TMath::Abs(pdgdau)==211){
 704         if(pdgdp<0 && pdgdau>0) return -1;
 705         if(pdgdp>0 && pdgdau<0) return -1;
 706         nPions++;
 707       }else if(TMath::Abs(pdgdau)==321){
 708         nKaons++;
 709       }else{
 710         return -1;
 711       }
 712     }
 713     if(nPions!=1) return -1;
 714     if(nKaons!=2) return -1;
 715     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 716       if(iDau<0) return -1;
 717     }
 718     if(isk0st) return 1;
 719     else if(isPhi) return 0;
 720     else return 3;
 721   }
 722   return -1;
 723 }
 724
 725 //______________________________________________________________________________
 726 Int_t AliAnalysisTaskCheckHFMCProd::CheckDstarDecay(Int_t labDstar, AliStack* stack) const{
 727
 728   if(labDstar<0) return -1;
 729   TParticle* dp = (TParticle*)stack->Particle(labDstar);
 730   Int_t pdgdp=dp->GetPdgCode();
 731   Int_t nDau=dp->GetNDaughters();
 732   if(nDau!=2) return -1;
 733
 734   Int_t nKaons=0;
 735   Int_t nPions=0;
 736   for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 737     if(iDau<0) return -1;
 738     TParticle* dau=(TParticle*)stack->Particle(iDau);
 739     Int_t pdgdau=dau->GetPdgCode();
 740     if(TMath::Abs(pdgdau)==421){
 741       for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 742         if(resDau<0) return -1;
 743         TParticle* resdau=(TParticle*)stack->Particle(resDau);
 744         Int_t pdgresdau=resdau->GetPdgCode();
 745         if(TMath::Abs(pdgresdau)==321){
 746           if(pdgdp>0 && pdgresdau>0) return -1;
 747           if(pdgdp<0 && pdgresdau<0) return -1;
 748           nKaons++;
 749         }
 750         if(TMath::Abs(pdgresdau)==211){
 751           if(pdgdp<0 && pdgresdau>0) return -1;
 752           if(pdgdp>0 && pdgresdau<0) return -1;
 753           nPions++;
 754         }
 755       }
 756     }else if(TMath::Abs(pdgdau)==211){
 757       if(pdgdp<0 && pdgdau>0) return -1;
 758       if(pdgdp>0 && pdgdau<0) return -1;
 759       nPions++;
 760     }else{
 761       return -1;
 762     }
 763   }
 764   if(nPions!=2) return -1;
 765   if(nKaons!=1) return -1;
 766   for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 767     if(iDau<0) return -1;
 768   }
 769   return 0;
 770
 771 }
 772
 773 //______________________________________________________________________________
 774 Int_t AliAnalysisTaskCheckHFMCProd::CheckLcDecay(Int_t labLc, AliStack* stack) const{
 775   if(labLc<0) return -1;
 776   TParticle* dp = (TParticle*)stack->Particle(labLc);
 777   Int_t pdgdp=dp->GetPdgCode();
 778   Int_t nDau=dp->GetNDaughters();
 779
 780   if(nDau==3){
 781     Int_t nKaons=0;
 782     Int_t nPions=0;
 783     Int_t nProtons=0;
 784     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 785       if(iDau<0) return -1;
 786       TParticle* dau=(TParticle*)stack->Particle(iDau);
 787       Int_t pdgdau=dau->GetPdgCode();
 788       if(TMath::Abs(pdgdau)==321){
 789         if(pdgdp>0 && pdgdau>0) return -1;
 790         if(pdgdp<0 && pdgdau<0) return -1;
 791         nKaons++;
 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 if(TMath::Abs(pdgdau)==2212){
 797         if(pdgdp<0 && pdgdau>0) return -1;
 798         if(pdgdp>0 && pdgdau<0) return -1;
 799         nProtons++;
 800       }else{
 801         return -1;
 802       }
 803     }
 804     if(nPions!=1) return -1;
 805     if(nKaons!=1) return -1;
 806     if(nProtons!=1) return -1;
 807     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 808       if(iDau<0) return -1;
 809     }
 810     return 0;
 811   }
 812
 813   if(nDau==2){
 814     Int_t nKaons=0;
 815     Int_t nPions=0;
 816     Int_t nProtons=0;
 817     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 818       if(iDau<0) return -1;
 819       TParticle* dau=(TParticle*)stack->Particle(iDau);
 820       Int_t pdgdau=dau->GetPdgCode();
 821       if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 || TMath::Abs(pdgdau)==2224
 822          || TMath::Abs(pdgdau)==3122 || TMath::Abs(pdgdau)==311){
 823         Int_t nDauRes=dau->GetNDaughters();
 824         if(nDauRes!=2)  return -1;
 825         for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){
 826           if(resDau<0) return -1;
 827           TParticle* resdau=(TParticle*)stack->Particle(resDau);
 828           Int_t pdgresdau=resdau->GetPdgCode();
 829           if(TMath::Abs(pdgresdau)==321){
 830             if(pdgdp>0 && pdgresdau>0) return -1;
 831             if(pdgdp<0 && pdgresdau<0) return -1;
 832             nKaons++;
 833           }
 834           if(TMath::Abs(pdgresdau)==211){
 835             if(pdgdp<0 && pdgresdau>0) return -1;
 836             if(pdgdp>0 && pdgresdau<0) return -1;
 837             nPions++;
 838           }
 839           if(TMath::Abs(pdgresdau)==2212){
 840             if(pdgdp<0 && pdgresdau>0) return -1;
 841             if(pdgdp>0 && pdgresdau<0) return -1;
 842             nProtons++;
 843           }
 844         }
 845       }else if(TMath::Abs(pdgdau)==321){
 846         if(pdgdp>0 && pdgdau>0) return -1;
 847         if(pdgdp<0 && pdgdau<0) return -1;
 848         nKaons++;
 849       }else if(TMath::Abs(pdgdau)==211){
 850         if(pdgdp<0 && pdgdau>0) return -1;
 851         if(pdgdp>0 && pdgdau<0) return -1;
 852         nPions++;
 853       }else if(TMath::Abs(pdgdau)==2212){
 854         if(pdgdp<0 && pdgdau>0) return -1;
 855         if(pdgdp>0 && pdgdau<0) return -1;
 856         nProtons++;
 857       }else{
 858         return -1;
 859       }
 860     }
 861     if(nPions!=1) return -1;
 862     if(nKaons!=1) return -1;
 863     if(nProtons!=1) return -1;
 864     for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
 865       if(iDau<0) return -1;
 866     }
 867     return 1;
 868   }
 869   return -1;
 870 }
 871