PWG3/vertexingHF/AliAnalysisTaskSELambdac.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /////////////////////////////////////////////////////////////
  17 //
  18 // AliAnalysisTaskSE for the extraction of signal(e.g Lambdac) of heavy flavor
  19 // decay candidates with the MC truth.
  20 // Authors: r.romita@gsi.de
  21 /////////////////////////////////////////////////////////////
  22
  23 #include <TClonesArray.h>
  24 #include <TNtuple.h>
  25 #include <TCanvas.h>
  26 #include <TList.h>
  27 #include <TString.h>
  28 #include <TH1F.h>
  29 #include <TH2F.h>
  30 #include <TDatabasePDG.h>
  31
  32 #include "AliAnalysisManager.h"
  33 #include "AliAODHandler.h"
  34 #include "AliAODEvent.h"
  35 #include "AliAODVertex.h"
  36 #include "AliAODTrack.h"
  37 #include "AliAODMCHeader.h"
  38 #include "AliAODMCParticle.h"
  39 #include "AliAODRecoDecayHF3Prong.h"
  40 #include "AliAnalysisVertexingHF.h"
  41 #include "AliAnalysisTaskSE.h"
  42 #include "AliAnalysisTaskSELambdac.h"
  43 #include "AliKFParticle.h"
  44 #include "AliAODPidHF.h"
  45 #include "AliRDHFCutsLctopKpi.h"
  46 #include "AliRDHFCuts.h"
  47 #include "AliKFVertex.h"
  48 #include "AliESDVertex.h"
  49
  50 ClassImp(AliAnalysisTaskSELambdac)
  51
  52
  53 //________________________________________________________________________
  54 AliAnalysisTaskSELambdac::AliAnalysisTaskSELambdac():
  55 AliAnalysisTaskSE(),
  56 fOutput(0),
  57 fHistNEvents(0),
  58 fhChi2(0),
  59 fhMassPtGreater3(0),
  60 fhMassPtGreater3TC(0),
  61 fNtupleLambdac(0),
  62 fUpmasslimit(2.486),
  63 fLowmasslimit(2.086),
  64 fNPtBins(0),
  65 fRDCutsAnalysis(0),
  66 fRDCutsProduction(0),
  67 fListCuts(0),
  68 fFillNtuple(kFALSE),
  69 fReadMC(kFALSE),
  70 fMCPid(kFALSE),
  71 fRealPid(kFALSE),
  72 fResPid(kTRUE),
  73 fUseKF(kFALSE),
  74 fVHF(0)
  75 {
  76    // Default constructor
  77 }
  78
  79 //________________________________________________________________________
  80 AliAnalysisTaskSELambdac::AliAnalysisTaskSELambdac(const char *name,Bool_t fillNtuple,AliRDHFCutsLctopKpi *lccutsana,AliRDHFCutsLctopKpi *lccutsprod):
  81 AliAnalysisTaskSE(name),
  82 fOutput(0),
  83 fHistNEvents(0),
  84 fhChi2(0),
  85 fhMassPtGreater3(0),
  86 fhMassPtGreater3TC(0),
  87 fNtupleLambdac(0),
  88 fUpmasslimit(2.486),
  89 fLowmasslimit(2.086),
  90 fNPtBins(0),
  91 fRDCutsAnalysis(lccutsana),
  92 fRDCutsProduction(lccutsprod),
  93 fListCuts(0),
  94 fFillNtuple(fillNtuple),
  95 fReadMC(kFALSE),
  96 fMCPid(kFALSE),
  97 fRealPid(kTRUE),
  98 fResPid(kFALSE),
  99 fUseKF(kFALSE),
 100 fVHF(0)
 101 {
 102    SetPtBinLimit(fRDCutsAnalysis->GetNPtBins()+1,fRDCutsAnalysis->GetPtBinLimits());
 103   // Default constructor
 104    // Output slot #1 writes into a TList container
 105   DefineOutput(1,TList::Class());  //My private output
 106   DefineOutput(2,TList::Class());
 107
 108   if(fFillNtuple){
 109     // Output slot #2 writes into a TNtuple container
 110     DefineOutput(3,TNtuple::Class());  //My private output
 111   }
 112 }
 113
 114 //________________________________________________________________________
 115 AliAnalysisTaskSELambdac::~AliAnalysisTaskSELambdac()
 116 {
 117   // Destructor
 118   if (fOutput) {
 119     delete fOutput;
 120     fOutput = 0;
 121   }
 122   if (fVHF) {
 123     delete fVHF;
 124     fVHF = 0;
 125   }
 126
 127  if(fRDCutsAnalysis){
 128     delete fRDCutsAnalysis;
 129     fRDCutsAnalysis = 0;
 130  }
 131  if(fRDCutsProduction){
 132     delete fRDCutsProduction;
 133     fRDCutsProduction = 0;
 134  }
 135
 136  if (fListCuts) {
 137    delete fListCuts;
 138    fListCuts = 0;
 139  }
 140 }
 141 //_________________________________________________________________
 142 void  AliAnalysisTaskSELambdac::SetMassLimits(Float_t range){
 143   fUpmasslimit = 2.286+range;
 144   fLowmasslimit = 2.286-range;
 145 }
 146 //_________________________________________________________________
 147 void  AliAnalysisTaskSELambdac::SetMassLimits(Float_t lowlimit, Float_t uplimit){
 148   if(uplimit>lowlimit)
 149     {
 150       fUpmasslimit = lowlimit;
 151       fLowmasslimit = uplimit;
 152     }
 153 }
 154
 155
 156 //________________________________________________________________________
 157 void AliAnalysisTaskSELambdac::SetPtBinLimit(Int_t n, Float_t* lim){
 158   // define pt bins for analysis
 159   if(n>kMaxPtBins){
 160     printf("Max. number of Pt bins = %d\n",kMaxPtBins);
 161     fNPtBins=kMaxPtBins;
 162     fArrayBinLimits[0]=0.;
 163     fArrayBinLimits[1]=2.;
 164     fArrayBinLimits[2]=3.;
 165     fArrayBinLimits[3]=4.;
 166     for(Int_t i=4; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
 167   }else{
 168     fNPtBins=n-1;
 169     fArrayBinLimits[0]=lim[0];
 170     for(Int_t i=1; i<fNPtBins+1; i++)
 171       if(lim[i]>fArrayBinLimits[i-1]){
 172         fArrayBinLimits[i]=lim[i];
 173       }
 174       else {
 175         fArrayBinLimits[i]=fArrayBinLimits[i-1];
 176       }
 177     for(Int_t i=fNPtBins; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
 178   }
 179   if(fDebug > 1){
 180     printf("Number of Pt bins = %d\n",fNPtBins);
 181     for(Int_t i=0; i<fNPtBins; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fArrayBinLimits[i],fArrayBinLimits[i+1]);
 182   }
 183 }
 184 //_________________________________________________________________
 185 Double_t  AliAnalysisTaskSELambdac::GetPtBinLimit(Int_t ibin) const{
 186   if(ibin>fNPtBins)return -1;
 187   return fArrayBinLimits[ibin];
 188 }
 189
 190 //_________________________________________________________________
 191 void AliAnalysisTaskSELambdac::Init()
 192 {
 193   // Initialization
 194
 195   if(fDebug > 1) printf("AnalysisTaskSELambdac::Init() \n");
 196
 197   fListCuts=new TList();
 198
 199   fListCuts->Add(new AliRDHFCutsLctopKpi(*fRDCutsAnalysis));
 200   fListCuts->Add(new AliRDHFCutsLctopKpi(*fRDCutsProduction));
 201   PostData(2,fListCuts);
 202   return;
 203 }
 204
 205 //________________________________________________________________________
 206 void AliAnalysisTaskSELambdac::UserCreateOutputObjects()
 207 {
 208   // Create the output container
 209   //
 210   if(fDebug > 1) printf("AnalysisTaskSELambdac::UserCreateOutputObjects() \n");
 211
 212   // Several histograms are more conveniently managed in a TList
 213   fOutput = new TList();
 214   fOutput->SetOwner();
 215   fOutput->SetName("OutputHistos");
 216
 217   TString hisname;
 218   Int_t index=0;
 219   Int_t indexLS=0;
 220   for(Int_t i=0;i<fNPtBins;i++){
 221
 222     index=GetHistoIndex(i);
 223     indexLS=GetLSHistoIndex(i);
 224
 225     hisname.Form("hMassPt%d",i);
 226     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 227     fMassHist[index]->Sumw2();
 228     hisname.Form("hCosPAllPt%d",i);
 229     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
 230     fCosPHist[index]->Sumw2();
 231     hisname.Form("hDLenAllPt%d",i);
 232     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
 233     fDLenHist[index]->Sumw2();
 234     hisname.Form("hSumd02AllPt%d",i);
 235     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
 236     fSumd02Hist[index]->Sumw2();
 237     hisname.Form("hSigVertAllPt%d",i);
 238     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 239     fSigVertHist[index]->Sumw2();
 240     hisname.Form("hPtMaxAllPt%d",i);
 241     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
 242     fPtMaxHist[index]->Sumw2();
 243
 244     hisname.Form("hDCAAllPt%d",i);
 245     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 246     fDCAHist[index]->Sumw2();
 247
 248
 249
 250     hisname.Form("hMassPt%dTC",i);
 251     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 252     fMassHistTC[index]->Sumw2();
 253
 254
 255
 256
 257
 258     hisname.Form("hCosPAllPt%dLS",i);
 259     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
 260     fCosPHistLS[index]->Sumw2();
 261     hisname.Form("hDLenAllPt%dLS",i);
 262     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
 263     fDLenHistLS[index]->Sumw2();
 264     hisname.Form("hSumd02AllPt%dLS",i);
 265     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
 266     fSumd02HistLS[index]->Sumw2();
 267     hisname.Form("hSigVertAllPt%dLS",i);
 268     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 269     fSigVertHistLS[index]->Sumw2();
 270     hisname.Form("hPtMaxAllPt%dLS",i);
 271     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
 272     fPtMaxHistLS[index]->Sumw2();
 273
 274     hisname.Form("hDCAAllPt%dLS",i);
 275     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 276     fDCAHistLS[index]->Sumw2();
 277
 278     hisname.Form("hLSPt%dLC",i);
 279     fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 280     fMassHistLS[indexLS]->Sumw2();
 281
 282     hisname.Form("hLSPt%dTC",i);
 283     fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 284     fMassHistLSTC[indexLS]->Sumw2();
 285
 286
 287
 288     index=GetSignalHistoIndex(i);
 289     indexLS++;
 290     hisname.Form("hSigPt%d",i);
 291     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 292     fMassHist[index]->Sumw2();
 293     hisname.Form("hCosPSigPt%d",i);
 294     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
 295     fCosPHist[index]->Sumw2();
 296     hisname.Form("hDLenSigPt%d",i);
 297     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
 298     fDLenHist[index]->Sumw2();
 299     hisname.Form("hSumd02SigPt%d",i);
 300     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
 301     fSumd02Hist[index]->Sumw2();
 302     hisname.Form("hSigVertSigPt%d",i);
 303     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 304     fSigVertHist[index]->Sumw2();
 305     hisname.Form("hPtMaxSigPt%d",i);
 306     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
 307     fPtMaxHist[index]->Sumw2();
 308
 309     hisname.Form("hDCASigPt%d",i);
 310     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 311     fDCAHist[index]->Sumw2();
 312
 313
 314     hisname.Form("hSigPt%dTC",i);
 315     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 316     fMassHistTC[index]->Sumw2();
 317
 318     hisname.Form("hLSPt%dLCnw",i);
 319     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 320     fMassHistLS[indexLS]->Sumw2();
 321     hisname.Form("hLSPt%dTCnw",i);
 322     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 323     fMassHistLSTC[indexLS]->Sumw2();
 324
 325
 326
 327     hisname.Form("hCosPSigPt%dLS",i);
 328     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
 329     fCosPHistLS[index]->Sumw2();
 330     hisname.Form("hDLenSigPt%dLS",i);
 331     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
 332     fDLenHistLS[index]->Sumw2();
 333     hisname.Form("hSumd02SigPt%dLS",i);
 334     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
 335     fSumd02HistLS[index]->Sumw2();
 336     hisname.Form("hSigVertSigPt%dLS",i);
 337     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 338     fSigVertHistLS[index]->Sumw2();
 339     hisname.Form("hPtMaxSigPt%dLS",i);
 340     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
 341     fPtMaxHistLS[index]->Sumw2();
 342
 343     hisname.Form("hDCASigPt%dLS",i);
 344     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 345     fDCAHistLS[index]->Sumw2();
 346
 347
 348
 349     index=GetBackgroundHistoIndex(i);
 350     indexLS++;
 351     hisname.Form("hBkgPt%d",i);
 352     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 353     fMassHist[index]->Sumw2();
 354     hisname.Form("hCosPBkgPt%d",i);
 355     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
 356     fCosPHist[index]->Sumw2();
 357     hisname.Form("hDLenBkgPt%d",i);
 358     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
 359     fDLenHist[index]->Sumw2();
 360     hisname.Form("hSumd02BkgPt%d",i);
 361     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
 362     fSumd02Hist[index]->Sumw2();
 363     hisname.Form("hSigVertBkgPt%d",i);
 364     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 365     fSigVertHist[index]->Sumw2();
 366     hisname.Form("hPtMaxBkgPt%d",i);
 367     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
 368     fPtMaxHist[index]->Sumw2();
 369
 370     hisname.Form("hDCABkgPt%d",i);
 371     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 372     fDCAHist[index]->Sumw2();
 373
 374
 375     hisname.Form("hBkgPt%dTC",i);
 376     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 377     fMassHistTC[index]->Sumw2();
 378
 379     hisname.Form("hLSPt%dLCntrip",i);
 380     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 381     fMassHistLS[indexLS]->Sumw2();
 382     hisname.Form("hLSPt%dTCntrip",i);
 383     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 384     fMassHistLSTC[indexLS]->Sumw2();
 385
 386
 387     hisname.Form("hCosPBkgPt%dLS",i);
 388     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
 389     fCosPHistLS[index]->Sumw2();
 390     hisname.Form("hDLenBkgPt%dLS",i);
 391     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
 392     fDLenHistLS[index]->Sumw2();
 393     hisname.Form("hSumd02BkgPt%dLS",i);
 394     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
 395     fSumd02HistLS[index]->Sumw2();
 396     hisname.Form("hSigVertBkgPt%dLS",i);
 397     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 398     fSigVertHistLS[index]->Sumw2();
 399     hisname.Form("hPtMaxBkgPt%dLS",i);
 400     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
 401     fPtMaxHistLS[index]->Sumw2();
 402     hisname.Form("hDCABkgPt%dLS",i);
 403     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 404     fDCAHistLS[index]->Sumw2();
 405
 406
 407     indexLS++;
 408     hisname.Form("hLSPt%dLCntripsinglecut",i);
 409     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 410     fMassHistLS[indexLS]->Sumw2();
 411     hisname.Form("hLSPt%dTCntripsinglecut",i);
 412     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 413     fMassHistLSTC[indexLS]->Sumw2();
 414
 415     indexLS++;
 416     hisname.Form("hLSPt%dLCspc",i);
 417     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 418     fMassHistLS[indexLS]->Sumw2();
 419     hisname.Form("hLSPt%dTCspc",i);
 420     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 421     fMassHistLSTC[indexLS]->Sumw2();
 422   }
 423
 424
 425   for(Int_t i=0; i<3*fNPtBins; i++){
 426     fOutput->Add(fMassHist[i]);
 427     fOutput->Add(fMassHistTC[i]);
 428     fOutput->Add(fCosPHist[i]);
 429     fOutput->Add(fDLenHist[i]);
 430     fOutput->Add(fSumd02Hist[i]);
 431     fOutput->Add(fSigVertHist[i]);
 432     fOutput->Add(fPtMaxHist[i]);
 433     fOutput->Add(fDCAHist[i]);
 434   }
 435
 436   fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
 437   fHistNEvents->Sumw2();
 438   fHistNEvents->SetMinimum(0);
 439   fOutput->Add(fHistNEvents);
 440
 441   fhChi2 = new TH1F("fhChi2", "Chi2",100,0.,10.);
 442   fhChi2->Sumw2();
 443   fOutput->Add(fhChi2);
 444
 445   fhMassPtGreater3=new TH1F("fhMassPtGreater3","Pt > 3 GeV/c",100,fLowmasslimit,fUpmasslimit);
 446   fhMassPtGreater3->Sumw2();
 447   fOutput->Add(fhMassPtGreater3);
 448   fhMassPtGreater3TC=new TH1F("fhMassPtGreater3TC","Pt > 3 GeV/c",100,fLowmasslimit,fUpmasslimit);
 449   fhMassPtGreater3TC->Sumw2();
 450   fOutput->Add(fhMassPtGreater3TC);
 451
 452
 453
 454   if(fFillNtuple){
 455     //OpenFile(3); // 2 is the slot number of the ntuple
 456
 457     fNtupleLambdac = new TNtuple("fNtupleLambdac","D +","pdg:Px:Py:Pz:PtTrue:VxTrue:VyTrue:VzTrue:Ptpi:PtK:Ptpi2:PtRec:PointingAngle:DecLeng:VxRec:VyRec:VzRec:InvMass:sigvert:d0Pi:d0K:d0Pi2:dca:d0square");
 458
 459   }
 460
 461   return;
 462 }
 463
 464 //________________________________________________________________________
 465 void AliAnalysisTaskSELambdac::UserExec(Option_t */*option*/)
 466 {
 467   // Execute analysis for current event:
 468   // heavy flavor candidates association to MC truth
 469
 470    AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 471    //tmp
 472  fHistNEvents->Fill(0); // count event
 473   // Post the data already here
 474   PostData(1,fOutput);
 475
 476   TClonesArray *array3Prong = 0;
 477   TClonesArray *arrayLikeSign =0;
 478   if(!aod && AODEvent() && IsStandardAOD()) {
 479     // In case there is an AOD handler writing a standard AOD, use the AOD
 480     // event in memory rather than the input (ESD) event.
 481     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 482     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 483     // have to taken from the AOD event hold by the AliAODExtension
 484     AliAODHandler* aodHandler = (AliAODHandler*)
 485       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 486     if(aodHandler->GetExtensions()) {
 487       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 488       AliAODEvent *aodFromExt = ext->GetAOD();
 489       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 490       arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
 491     }
 492   } else if(aod) {
 493     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 494     arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
 495   }
 496
 497   if(!array3Prong || !aod) {
 498     printf("AliAnalysisTaskSELambdac::UserExec: Charm3Prong branch not found!\n");
 499     return;
 500   }
 501   if(!arrayLikeSign) {
 502     printf("AliAnalysisTaskSELambdac::UserExec: LikeSign3Prong branch not found!\n");
 503   //  return;
 504   }
 505
 506
 507   TClonesArray *arrayMC=0;
 508   AliAODMCHeader *mcHeader=0;
 509
 510   // AOD primary vertex
 511   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 512   if(!vtx1) return;
 513   if(vtx1==0x0) return;
 514
 515   if(!fRDCutsProduction->IsEventSelected(aod)) return;
 516
 517   // load MC particles
 518   if(fReadMC){
 519
 520     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 521     if(!arrayMC) {
 522       printf("AliAnalysisTaskSELambdac::UserExec: MC particles branch not found!\n");
 523       return;
 524     }
 525
 526   // load MC header
 527     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 528     if(!mcHeader) {
 529     printf("AliAnalysisTaskSELambdac::UserExec: MC header branch not found!\n");
 530     return;
 531     }
 532   }
 533
 534   Int_t n3Prong = array3Prong->GetEntriesFast();
 535
 536
 537   Int_t nOS=0;
 538   Int_t index;
 539   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 540     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 541
 542
 543     Bool_t unsetvtx=kFALSE;
 544     if(!d->GetOwnPrimaryVtx()){
 545       d->SetOwnPrimaryVtx(vtx1);
 546       unsetvtx=kTRUE;
 547     }
 548
 549     Int_t isSelectedTracks = fRDCutsProduction->IsSelected(d,AliRDHFCuts::kTracks);
 550     if(!isSelectedTracks) continue;
 551
 552     Int_t selection=fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod);
 553     if(selection>0) {
 554       Int_t iPtBin = -1;
 555       Double_t ptCand = d->Pt();
 556
 557       for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
 558         if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
 559       }
 560
 561       Int_t labDp=-1;
 562       Float_t deltaPx=0.;
 563       Float_t deltaPy=0.;
 564       Float_t deltaPz=0.;
 565       Float_t truePt=0.;
 566       Float_t xDecay=0.;
 567       Float_t yDecay=0.;
 568       Float_t zDecay=0.;
 569       Float_t pdgCode=-2;
 570       Bool_t isSignal=kFALSE;
 571       Float_t pdgCode1=-1;
 572       Float_t pdgCode2=-1;
 573       if(fReadMC){
 574         labDp = MatchToMCLambdac(d,arrayMC);
 575         //
 576         if(labDp>0){
 577         isSignal=kTRUE;
 578           AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
 579           AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
 580           AliAODMCParticle *dg1 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(1));
 581           deltaPx=partDp->Px()-d->Px();
 582           deltaPy=partDp->Py()-d->Py();
 583           deltaPz=partDp->Pz()-d->Pz();
 584           truePt=partDp->Pt();
 585           xDecay=dg0->Xv();
 586           yDecay=dg0->Yv();
 587           zDecay=dg0->Zv();
 588           pdgCode=TMath::Abs(partDp->GetPdgCode());
 589           pdgCode1=TMath::Abs(dg0->GetPdgCode());
 590           pdgCode2=TMath::Abs(dg1->GetPdgCode());
 591         }else{
 592           pdgCode=-1;
 593         }
 594       }
 595
 596       Double_t invMasspKpi=-1.;
 597       Double_t invMasspiKp=-1.;
 598       Int_t pdgs[3]={0,0,0};
 599       Double_t field=aod->GetMagneticField();
 600       //apply MC PID
 601       if(fReadMC && fMCPid){
 602
 603        if(IspKpiMC(d,arrayMC)) {
 604         invMasspKpi=d->InvMassLcpKpi();
 605         if(fUseKF){
 606          pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
 607          if(!VertexingKF(d,pdgs,field)) invMasspKpi=-1.;
 608         }
 609        }
 610        if(IspiKpMC(d,arrayMC)) {
 611         invMasspiKp=d->InvMassLcpiKp();
 612         if(fUseKF){
 613          pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
 614          if(!VertexingKF(d,pdgs,field)) invMasspiKp=-1.;
 615         }
 616        }
 617       }
 618       // apply realistic PID
 619       if(fRealPid){
 620        if(selection==1 || selection==3) {
 621         invMasspKpi=d->InvMassLcpKpi();
 622         pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
 623         if(fUseKF){
 624          pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
 625          if(!VertexingKF(d,pdgs,field)) invMasspKpi=-1.;
 626         }
 627        }
 628        if(selection>=2) {
 629         invMasspiKp=d->InvMassLcpiKp();
 630          pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
 631         if(fUseKF){
 632          pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
 633          if(!VertexingKF(d,pdgs,field)) invMasspiKp=-1.;
 634         }
 635
 636        }
 637       }
 638
 639       //apply PID using resonances
 640       if(fResPid){
 641        if(IspKpiResonant(d,field) && (selection==3 || selection==1)) {
 642         invMasspKpi=d->InvMassLcpKpi();
 643         if(fUseKF){
 644          pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
 645          if(!VertexingKF(d,pdgs,field)) invMasspKpi=-1.;
 646         }
 647        }
 648        if(IspiKpResonant(d,field) && selection>=2) {
 649         invMasspiKp=d->InvMassLcpiKp();
 650         if(fUseKF){
 651          pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
 652          if(!VertexingKF(d,pdgs,field)) invMasspiKp=-1.;
 653         }
 654        }
 655       }
 656       // no PID
 657       if(!fResPid && !fRealPid && !fMCPid){
 658        if(selection==2 || selection==3) invMasspiKp=d->InvMassLcpiKp();
 659        if(fUseKF){
 660          pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
 661          if(!VertexingKF(d,pdgs,field)) invMasspiKp=-1.;
 662         }
 663        if(selection==1 || selection==3) invMasspKpi=d->InvMassLcpKpi();
 664        if(fUseKF){
 665         pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
 666         if(!VertexingKF(d,pdgs,field)) invMasspKpi=-1.;
 667        }
 668       }
 669
 670       Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
 671
 672       if(invMasspiKp<0. && invMasspKpi<0.) continue;
 673
 674
 675       Float_t tmp[24];
 676       if(fFillNtuple){
 677         tmp[0]=pdgCode;
 678         tmp[1]=deltaPx;
 679         tmp[2]=deltaPy;
 680         tmp[3]=deltaPz;
 681         tmp[4]=truePt;
 682         tmp[5]=xDecay;
 683         tmp[6]=yDecay;
 684         tmp[7]=zDecay;
 685         if(pdgCode1==2212) {tmp[8]=d->PtProng(0);}else{tmp[8]=0.;}
 686         if(pdgCode1==211) {tmp[10]=d->PtProng(0);}else{tmp[10]=0.;}
 687         tmp[9]=d->PtProng(1);
 688         if(pdgCode2==211) {tmp[10]=d->PtProng(2);}else{tmp[10]=0.;}
 689         tmp[11]=d->Pt();
 690         tmp[12]=d->CosPointingAngle();
 691         tmp[13]=d->DecayLength();
 692         tmp[14]=d->Xv();
 693         tmp[15]=d->Yv();
 694         tmp[16]=d->Zv();
 695         if(invMasspiKp>0.) tmp[17]=invMasspiKp;
 696         if(invMasspKpi>0.) tmp[17]=invMasspKpi;
 697         tmp[18]=d->GetSigmaVert();
 698         tmp[19]=d->Getd0Prong(0);
 699         tmp[20]=d->Getd0Prong(1);
 700         tmp[21]=d->Getd0Prong(2);
 701         tmp[22]=d->GetDCA();
 702         tmp[23]=d->Prodd0d0();
 703         fNtupleLambdac->Fill(tmp);
 704         PostData(3,fNtupleLambdac);
 705       }
 706       Double_t dlen=d->DecayLength();
 707       Double_t cosp=d->CosPointingAngle();
 708       Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 709       Double_t dca=d->GetDCA();
 710 Double_t ptmax=0;
 711       for(Int_t i=0;i<3;i++){
 712         if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
 713       }
 714      if(d->Pt()>3.){
 715         if(invMasspiKp>0. && invMasspKpi>0.){
 716         if(invMasspiKp>0.) fhMassPtGreater3->Fill(invMasspiKp,0.5);
 717         if(invMasspKpi>0.) fhMassPtGreater3->Fill(invMasspKpi,0.5);
 718         }else{
 719          if(invMasspiKp>0.) fhMassPtGreater3->Fill(invMasspiKp);
 720          if(invMasspKpi>0.) fhMassPtGreater3->Fill(invMasspKpi);
 721         }
 722         if(passTightCuts>0){
 723          if(invMasspiKp>0. && invMasspKpi>0.){
 724          if(invMasspiKp>0.) fhMassPtGreater3TC->Fill(invMasspiKp,0.5);
 725          if(invMasspKpi>0.) fhMassPtGreater3TC->Fill(invMasspKpi,0.5);
 726          }else{
 727           if(invMasspiKp>0.) fhMassPtGreater3TC->Fill(invMasspiKp);
 728           if(invMasspKpi>0.) fhMassPtGreater3TC->Fill(invMasspKpi);
 729          }
 730         }
 731       }
 732       if(iPtBin>=0){
 733
 734         index=GetHistoIndex(iPtBin);
 735         if(invMasspiKp>0. && invMasspKpi>0.){
 736         if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp,0.5);
 737         if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi,0.5);
 738         }else{
 739          if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
 740          if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
 741         }
 742
 743         fCosPHist[index]->Fill(cosp);
 744         fDLenHist[index]->Fill(dlen);
 745         fSumd02Hist[index]->Fill(sumD02);
 746         fPtMaxHist[index]->Fill(ptmax);
 747         fDCAHist[index]->Fill(dca);
 748
 749         if(passTightCuts>0){
 750          if(invMasspiKp>0. && invMasspKpi>0. && passTightCuts==3){
 751           if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp,0.5);
 752           if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi,0.5);
 753          }else{
 754           if(invMasspiKp>0. && passTightCuts==2) fMassHistTC[index]->Fill(invMasspiKp);
 755           if(invMasspKpi>0. && passTightCuts==1) fMassHistTC[index]->Fill(invMasspKpi);
 756          }
 757         }
 758
 759         if(fReadMC){
 760           if(labDp>0) {
 761             index=GetSignalHistoIndex(iPtBin);
 762             if(invMasspiKp>0. && invMasspKpi>0.){
 763              if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp,0.5);
 764              if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi,0.5);
 765             }else{
 766              if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
 767              if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
 768             }
 769             fCosPHist[index]->Fill(cosp);
 770             fDLenHist[index]->Fill(dlen);
 771             fSumd02Hist[index]->Fill(sumD02);
 772             fPtMaxHist[index]->Fill(ptmax);
 773             fDCAHist[index]->Fill(dca);
 774             if(passTightCuts>0){
 775              if(invMasspiKp>0. && invMasspKpi>0. && passTightCuts==3){
 776               if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp,0.5);
 777               if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi,0.5);
 778              }else{
 779               if(invMasspiKp>0. && passTightCuts==2) fMassHistTC[index]->Fill(invMasspiKp);
 780               if(invMasspKpi>0.&& passTightCuts==1) fMassHistTC[index]->Fill(invMasspKpi);
 781              }
 782             }
 783
 784           }else{
 785             index=GetBackgroundHistoIndex(iPtBin);
 786             if(invMasspiKp>0. && invMasspKpi>0.){
 787              fMassHist[index]->Fill(invMasspiKp,0.5);
 788              fMassHist[index]->Fill(invMasspKpi,0.5);
 789             }else{
 790              if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
 791              if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
 792            }
 793             fCosPHist[index]->Fill(cosp);
 794             fDLenHist[index]->Fill(dlen);
 795             fSumd02Hist[index]->Fill(sumD02);
 796             fPtMaxHist[index]->Fill(ptmax);
 797             fDCAHist[index]->Fill(dca);
 798             if(invMasspiKp>0. && invMasspKpi>0. && passTightCuts==3){
 799                fMassHistTC[index]->Fill(invMasspiKp,0.5);
 800                fMassHistTC[index]->Fill(invMasspKpi,0.5);
 801               }else{
 802                if(invMasspiKp>0. && passTightCuts==2) fMassHistTC[index]->Fill(invMasspiKp);
 803                if(invMasspKpi>0. && passTightCuts==1) fMassHistTC[index]->Fill(invMasspKpi);
 804               }
 805           }
 806         }
 807       }
 808       /*
 809       //start OS analysis
 810       if(labDp<0)fHistOSbkg->Fill(d->InvMassDplus());
 811       fHistOS->Fill(d->InvMassDplus());
 812       */
 813       nOS++;
 814     }
 815     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 816   }
 817
 818   PostData(1,fOutput);
 819   return;
 820 }
 821
 822
 823
 824 //________________________________________________________________________
 825 void AliAnalysisTaskSELambdac::Terminate(Option_t */*option*/)
 826 {
 827   // Terminate analysis
 828   //
 829   if(fDebug > 1) printf("AnalysisTaskSELambdac: Terminate() \n");
 830
 831   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 832   if (!fOutput) {
 833     printf("ERROR: fOutput not available\n");
 834     return;
 835   }
 836  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
 837
 838  TString hisname;
 839  Int_t index=0;
 840
 841
 842  for(Int_t i=0;i<fNPtBins;i++){
 843     index=GetHistoIndex(i);
 844     hisname.Form("hMassPt%d",i);
 845     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 846      hisname.Form("hCosPAllPt%d",i);
 847     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 848      hisname.Form("hDLenAllPt%d",i);
 849     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 850      hisname.Form("hSumd02AllPt%d",i);
 851      fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 852      hisname.Form("hSigVertAllPt%d",i);
 853      fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 854      hisname.Form("hPtMaxAllPt%d",i);
 855      fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 856      hisname.Form("hDCAAllPt%d",i);
 857      fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 858      hisname.Form("hMassPt%dTC",i);
 859     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 860
 861     index=GetSignalHistoIndex(i);
 862     hisname.Form("hSigPt%d",i);
 863     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 864     hisname.Form("hCosPSigPt%d",i);
 865     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 866     hisname.Form("hDLenSigPt%d",i);
 867     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 868     hisname.Form("hSumd02SigPt%d",i);
 869     fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 870     hisname.Form("hSigVertSigPt%d",i);
 871     fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 872     hisname.Form("hPtMaxSigPt%d",i);
 873     fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 874     hisname.Form("hDCASigPt%d",i);
 875     fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 876
 877     hisname.Form("hSigPt%dTC",i);
 878     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 879
 880     index=GetBackgroundHistoIndex(i);
 881     hisname.Form("hBkgPt%d",i);
 882     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 883     hisname.Form("hCosPBkgPt%d",i);
 884     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 885     hisname.Form("hDLenBkgPt%d",i);
 886     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 887     hisname.Form("hSumd02BkgPt%d",i);
 888     fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 889     hisname.Form("hSigVertBkgPt%d",i);
 890     fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 891     hisname.Form("hPtMaxBkgPt%d",i);
 892     fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 893     hisname.Form("hDCABkgPt%d",i);
 894     fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 895     hisname.Form("hBkgPt%dTC",i);
 896     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 897
 898  }
 899
 900   if(fFillNtuple){
 901     fNtupleLambdac = dynamic_cast<TNtuple*>(GetOutputData(3));
 902   }
 903
 904
 905  return;
 906 }
 907
 908 //________________________________________________________________________
 909 Int_t AliAnalysisTaskSELambdac::MatchToMCLambdac(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC) const{
 910   // check if the candidate is a Lambdac decaying in pKpi or in the resonant channels
 911   Int_t lambdacLab[3]={0,0,0};
 912   Int_t pdgs[3]={0,0,0};
 913   for(Int_t i=0;i<3;i++){
 914    AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
 915    Int_t lab=daugh->GetLabel();
 916    if(lab<0) return 0;
 917    AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab);
 918    if(!part) continue;
 919    pdgs[i]=part->GetPdgCode();
 920    Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
 921    if(partPdgcode==211 || partPdgcode==321 || partPdgcode==2212){
 922         Int_t motherLabel=part->GetMother();
 923              if(motherLabel<0) return 0;
 924    AliAODMCParticle *motherPart = (AliAODMCParticle*)arrayMC->At(motherLabel);
 925    if(!motherPart) continue;
 926         Int_t motherPdg = TMath::Abs(motherPart->GetPdgCode());
 927    if(motherPdg==4122) {
 928         if(GetLambdacDaugh(motherPart,arrayMC)){lambdacLab[i]=motherLabel;continue;}
 929               }
 930    if(motherPdg==313 || motherPdg==2224 || motherPdg==3124){
 931          Int_t GmotherLabel=motherPart->GetMother();
 932                if(GmotherLabel<0) return 0;
 933          AliAODMCParticle *GmotherPart = (AliAODMCParticle*)arrayMC->At(GmotherLabel);
 934          if(!GmotherPart) continue;
 935                Int_t GmotherPdg = TMath::Abs(GmotherPart->GetPdgCode());
 936          if(GmotherPdg==4122) {
 937                 if(GetLambdacDaugh(GmotherPart,arrayMC)) {lambdacLab[i]=GmotherLabel;continue;}
 938               }
 939         }
 940      }
 941   }
 942
 943  if(lambdacLab[0]==lambdacLab[1] && lambdacLab[1]==lambdacLab[2]) {return lambdacLab[0];}
 944  return 0;
 945
 946 }
 947 //------------------------
 948 Bool_t AliAnalysisTaskSELambdac::GetLambdacDaugh(AliAODMCParticle *part,TClonesArray *arrayMC) const{
 949  // check if the particle is a lambdac and if its decay mode is the correct one
 950  Int_t numberOfLambdac=0;
 951  if(TMath::Abs(part->GetPdgCode())!=4122) return kFALSE;
 952  Int_t daugh_tmp[2];
 953  daugh_tmp[0]=part->GetDaughter(0);
 954  daugh_tmp[1]=part->GetDaughter(1);
 955  Int_t nDaugh = (Int_t)part->GetNDaughters();
 956  if(nDaugh<2) return kFALSE;
 957  if(nDaugh>3) return kFALSE;
 958  AliAODMCParticle* pdaugh1 = (AliAODMCParticle*)arrayMC->At(part->GetDaughter(0));
 959    if(!pdaugh1) {return kFALSE;}
 960  Int_t number1 = TMath::Abs(pdaugh1->GetPdgCode());
 961  AliAODMCParticle* pdaugh2 = (AliAODMCParticle*)arrayMC->At(part->GetDaughter(1));
 962    if(!pdaugh2) {return kFALSE;}
 963   Int_t number2 = TMath::Abs(pdaugh2->GetPdgCode());
 964
 965   if(nDaugh==3){
 966    Int_t thirdDaugh=part->GetDaughter(1)-1;
 967             AliAODMCParticle* pdaugh3 = (AliAODMCParticle*)arrayMC->At(thirdDaugh);
 968   Int_t number3 = TMath::Abs(pdaugh3->GetPdgCode());
 969    if((number1==321 && number2==211 && number3==2212) || (number1==211 && number2==321 && number3==2212) || (number1==211 && number2==2212 && number3==321) || (number1==321 && number2==2212 && number3==211) || (number1==2212 && number2==321 && number3==211) || (number1==2212 && number2==211 && number3==321)) {
 970    numberOfLambdac++;
 971    }
 972   }
 973
 974  if(nDaugh==2){
 975
 976   //Lambda resonant
 977
 978   //Lambda -> p K*0
 979   //
 980   Int_t nfiglieK=0;
 981
 982    if((number1==2212 && number2==313)){
 983      nfiglieK=pdaugh2->GetNDaughters();
 984      if(nfiglieK!=2) return kFALSE;
 985      AliAODMCParticle* pdaughK1 = (AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
 986      AliAODMCParticle* pdaughK2 = (AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
 987      if(!pdaughK1) return kFALSE;
 988      if(!pdaughK2) return kFALSE;
 989      if((TMath::Abs(pdaughK1->GetPdgCode())==211 && TMath::Abs(pdaughK2->GetPdgCode())==321) || (TMath::Abs(pdaughK1->GetPdgCode())==321 && TMath::Abs(pdaughK2->GetPdgCode())==211)) numberOfLambdac++;
 990     }
 991
 992    if((number1==313 && number2==2212)){
 993     nfiglieK=pdaugh1->GetNDaughters();
 994     if(nfiglieK!=2) return kFALSE;
 995     AliAODMCParticle* pdaughK1 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
 996     AliAODMCParticle* pdaughK2 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
 997     if(!pdaughK1) return kFALSE;
 998     if(!pdaughK2) return kFALSE;
 999      if((TMath::Abs(pdaughK1->GetPdgCode())==211 && TMath::Abs(pdaughK2->GetPdgCode())==321) || (TMath::Abs(pdaughK1->GetPdgCode())==321 && TMath::Abs(pdaughK2->GetPdgCode())==211)) numberOfLambdac++;
1000    }
1001
1002    //Lambda -> Delta++ k
1003    Int_t nfiglieDelta=0;
1004    if(number1==321 && number2==2224){
1005     nfiglieDelta=pdaugh2->GetNDaughters();
1006     if(nfiglieDelta!=2) return kFALSE;
1007     AliAODMCParticle *pdaughD1=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
1008     AliAODMCParticle *pdaughD2=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
1009     if(!pdaughD1) return kFALSE;
1010     if(!pdaughD2) return kFALSE;
1011     if((TMath::Abs(pdaughD1->GetPdgCode())==211 && TMath::Abs(pdaughD2->GetPdgCode())==2212) || (TMath::Abs(pdaughD1->GetPdgCode())==2212 && TMath::Abs(pdaughD2->GetPdgCode())==211)) numberOfLambdac++;
1012    }
1013    if(number1==2224 && number2==321){
1014     nfiglieDelta=pdaugh1->GetNDaughters();
1015     if(nfiglieDelta!=2) return kFALSE;
1016     AliAODMCParticle* pdaughD1 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
1017     AliAODMCParticle* pdaughD2 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
1018     if(!pdaughD1) return kFALSE;
1019     if(!pdaughD2) return kFALSE;
1020     if((TMath::Abs(pdaughD1->GetPdgCode())==211 && TMath::Abs(pdaughD2->GetPdgCode())==2212) || (TMath::Abs(pdaughD1->GetPdgCode())==2212 && TMath::Abs(pdaughD2->GetPdgCode())==211)) numberOfLambdac++;
1021    }
1022
1023
1024    //Lambdac -> Lambda(1520) pi
1025    Int_t nfiglieLa=0;
1026    if(number1==3124 && number2==211){
1027     nfiglieLa=pdaugh1->GetNDaughters();
1028     if(nfiglieLa!=2) return kFALSE;
1029     AliAODMCParticle *pdaughL1=(AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
1030     AliAODMCParticle *pdaughL2=(AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
1031     if(!pdaughL1) return kFALSE;
1032     if(!pdaughL2) return kFALSE;
1033     if((TMath::Abs(pdaughL1->GetPdgCode())==321 && TMath::Abs(pdaughL2->GetPdgCode())==2212) || (TMath::Abs(pdaughL1->GetPdgCode())==2212 && TMath::Abs(pdaughL2->GetPdgCode())==321)) numberOfLambdac++;
1034    }
1035    if(number1==211 && number2==3124){
1036     nfiglieLa=pdaugh2->GetNDaughters();
1037     if(nfiglieLa!=2) return kFALSE;
1038     AliAODMCParticle *pdaughL1=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
1039     AliAODMCParticle *pdaughL2=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
1040     if(!pdaughL1) return kFALSE;
1041     if(!pdaughL2) return kFALSE;
1042     if((TMath::Abs(pdaughL1->GetPdgCode())==321 && TMath::Abs(pdaughL2->GetPdgCode())==2212) || (TMath::Abs(pdaughL1->GetPdgCode())==2212 && TMath::Abs(pdaughL2->GetPdgCode())==321)) numberOfLambdac++;
1043
1044    }
1045  }
1046
1047  if(numberOfLambdac>0) {return kTRUE;}
1048          return kFALSE;
1049 }
1050 //-----------------------------
1051 Bool_t AliAnalysisTaskSELambdac::IspKpiMC(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC) const{
1052   // Apply MC PID
1053    Int_t lab[3]={0,0,0},pdgs[3]={0,0,0};
1054    for(Int_t i=0;i<3;i++){
1055     AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
1056     lab[i]=daugh->GetLabel();
1057     if(lab[i]<0) return kFALSE;
1058     AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab[i]);
1059     if(!part) return kFALSE;
1060     pdgs[i]=TMath::Abs(part->GetPdgCode());
1061    }
1062
1063    if(pdgs[0]==2212 && pdgs[1]==321 && pdgs[2]==211) return kTRUE;
1064
1065    return kFALSE;
1066 }
1067 //-----------------------------
1068 Bool_t AliAnalysisTaskSELambdac::IspiKpMC(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC) const{
1069
1070   // Apply MC PID
1071    Int_t lab[3]={0,0,0},pdgs[3]={0,0,0};
1072    for(Int_t i=0;i<3;i++){
1073     AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
1074     lab[i]=daugh->GetLabel();
1075     if(lab[i]<0) return kFALSE;
1076     AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab[i]);
1077     if(!part) return kFALSE;
1078     pdgs[i]=TMath::Abs(part->GetPdgCode());
1079    }
1080
1081    if(pdgs[2]==2212 && pdgs[1]==321 && pdgs[0]==211) {return kTRUE;}
1082
1083    return kFALSE;
1084 }
1085 //--------------------------------------
1086 Bool_t AliAnalysisTaskSELambdac::VertexingKF(AliAODRecoDecayHF3Prong *d,Int_t *pdgs,Double_t field) const{
1087  // apply vertexing KF
1088    Int_t iprongs[3]={0,1,2};
1089    Double_t mass[2]={0.,0.};
1090  //topological constr
1091    AliKFParticle *lambdac=d->ApplyVertexingKF(iprongs,3,pdgs,kTRUE,field,mass);
1092    if(!lambdac) return kFALSE;
1093 //  Double_t probTot=TMath::Prob(Lambdac->GetChi2(),Lambdac->GetNDF());
1094 //  if(probTot<fCutsKF[0]) return kFALSE;
1095   if(lambdac->GetChi2()>fCutsKF[0]) return kFALSE;
1096  //mass constr for K*
1097    Int_t ipRes[2];
1098    Int_t pdgres[2];
1099    mass[0]=0.8961;mass[1]=0.03;
1100    if(TMath::Abs(pdgs[0])==211){
1101     ipRes[0]=0;ipRes[1]=1;
1102     pdgres[0]=pdgs[0];pdgres[1]=321;
1103    }
1104    if(TMath::Abs(pdgs[2])==211){
1105     ipRes[0]=2;ipRes[1]=1;
1106     pdgres[0]=pdgs[2];pdgres[1]=321;
1107    }
1108    AliKFParticle *kappaStar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1109
1110   Double_t probKstar=TMath::Prob(kappaStar->GetChi2(),kappaStar->GetNDF());
1111   if(probKstar>fCutsKF[1]) {
1112     AliAODTrack *esdProng1=(AliAODTrack*)d->GetDaughter(ipRes[0]);
1113     AliAODTrack *esdProng2=(AliAODTrack*)d->GetDaughter(ipRes[1]);
1114     AliKFParticle prong1(*esdProng1,pdgres[0]);
1115     AliKFParticle prong2(*esdProng2,pdgres[1]);
1116     if(kappaStar->GetPt()<fCutsKF[2] && prong1.GetAngle(prong2)>fCutsKF[3]) return kFALSE;
1117   }
1118  //mass constr for Lambda
1119    mass[0]=1.520;mass[1]=0.005;
1120    if(TMath::Abs(pdgs[0])==2212){
1121     ipRes[0]=0;ipRes[1]=1;
1122     pdgres[0]=pdgs[0];pdgres[1]=pdgs[1];
1123    }
1124    if(TMath::Abs(pdgs[2])==2212){
1125     ipRes[0]=2;ipRes[1]=1;
1126     pdgres[0]=pdgs[2];pdgres[1]=pdgs[1];
1127    }
1128    AliKFParticle *lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1129   Double_t probLa=TMath::Prob(lambda1520->GetChi2(),lambda1520->GetNDF());
1130   if(probLa>fCutsKF[4]) {
1131     AliAODTrack *esdProng1=(AliAODTrack*)d->GetDaughter(ipRes[0]);
1132     AliAODTrack *esdProng2=(AliAODTrack*)d->GetDaughter(ipRes[1]);
1133     AliKFParticle prong1(*esdProng1,pdgres[0]);
1134     AliKFParticle prong2(*esdProng2,pdgres[1]);
1135     if(lambda1520->GetPt()<fCutsKF[5] && prong1.GetAngle(prong2)>fCutsKF[6]) return kFALSE;
1136   }
1137  //mass constr for Delta
1138    mass[0]=1.2;mass[1]=0.15;
1139    ipRes[0]=0;ipRes[1]=2;
1140    pdgres[0]=pdgs[0];pdgres[1]=pdgs[2];
1141    AliKFParticle *delta=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1142   Double_t probDelta=TMath::Prob(delta->GetChi2(),delta->GetNDF());
1143   if(probDelta>fCutsKF[7]) {
1144     AliAODTrack *esdProng1=(AliAODTrack*)d->GetDaughter(ipRes[0]);
1145     AliAODTrack *esdProng2=(AliAODTrack*)d->GetDaughter(ipRes[1]);
1146     AliKFParticle prong1(*esdProng1,pdgres[0]);
1147     AliKFParticle prong2(*esdProng2,pdgres[1]);
1148     if(delta->GetPt()<fCutsKF[8] && prong1.GetAngle(prong2)>fCutsKF[9]) return kFALSE;
1149   }
1150  return kTRUE;
1151 }
1152 //-------------------------------------
1153 Bool_t AliAnalysisTaskSELambdac::IspiKpResonant(AliAODRecoDecayHF3Prong *d,Double_t field) const{
1154
1155  // apply PID using the resonant channels
1156  //if lambda* -> pk
1157         Double_t mass[2]={1.520,0.005};
1158         Int_t ipRes[2]={1,2};
1159         Int_t pdgres[2]={321,2212};
1160         AliKFParticle *lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1161         Double_t probLa=TMath::Prob(lambda1520->GetChi2(),lambda1520->GetNDF());
1162         if(probLa>0.9) return kTRUE;
1163  //K* -> kpi
1164         mass[0]=0.8961;mass[1]=0.03;
1165         ipRes[0]=0;ipRes[1]=1;
1166         pdgres[0]=211;pdgres[1]=321;
1167         AliKFParticle *kstar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1168         Double_t probKa=TMath::Prob(kstar->GetChi2(),kstar->GetNDF());
1169        if(probKa>0.9) return kTRUE;
1170
1171  return kFALSE;
1172
1173 }
1174 //-------------------------------------
1175 Bool_t AliAnalysisTaskSELambdac::IspKpiResonant(AliAODRecoDecayHF3Prong *d,Double_t field) const{
1176
1177  // apply PID using the resonant channels
1178  //if lambda* -> pk
1179         Double_t mass[2]={1.520,0.005};
1180         Int_t ipRes[2]={0,1};
1181         Int_t pdgres[2]={2212,321};
1182         AliKFParticle *lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1183         Double_t probLa=TMath::Prob(lambda1520->GetChi2(),lambda1520->GetNDF());
1184         if(probLa>0.9) return kTRUE;
1185  //K* -> kpi
1186         mass[0]=0.8961;mass[1]=0.03;
1187         ipRes[0]=1;ipRes[1]=2;
1188         pdgres[1]=211;pdgres[0]=321;
1189         AliKFParticle *kstar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1190         Double_t probKa=TMath::Prob(kstar->GetChi2(),kstar->GetNDF());
1191         if(probKa>0.9) return kTRUE;
1192
1193  return kFALSE;
1194
1195 }
1196