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