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