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