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