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
 458   TClonesArray *array3Prong = 0;
 459   TClonesArray *arrayLikeSign =0;
 460   if(!aod && AODEvent() && IsStandardAOD()) {
 461     // In case there is an AOD handler writing a standard AOD, use the AOD
 462     // event in memory rather than the input (ESD) event.
 463     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 464     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 465     // have to taken from the AOD event hold by the AliAODExtension
 466     AliAODHandler* aodHandler = (AliAODHandler*)
 467       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 468     if(aodHandler->GetExtensions()) {
 469       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 470       AliAODEvent *aodFromExt = ext->GetAOD();
 471       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 472       arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
 473     }
 474   } else {
 475     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 476     arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
 477   }
 478
 479   if(!array3Prong) {
 480     printf("AliAnalysisTaskSELambdac::UserExec: Charm3Prong branch not found!\n");
 481     return;
 482   }
 483   if(!arrayLikeSign) {
 484     printf("AliAnalysisTaskSELambdac::UserExec: LikeSign3Prong branch not found!\n");
 485   //  return;
 486   }
 487
 488   // fix for temporary bug in ESDfilter
 489   // the AODs with null vertex pointer didn't pass the PhysSel
 490   if(!aod->GetPrimaryVertex()) return;
 491
 492   fHistNEvents->Fill(0); // count event
 493   // Post the data already here
 494   PostData(1,fOutput);
 495
 496   TClonesArray *arrayMC=0;
 497   AliAODMCHeader *mcHeader=0;
 498
 499   // AOD primary vertex
 500   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 501
 502   // load MC particles
 503   if(fReadMC){
 504
 505     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 506     if(!arrayMC) {
 507       printf("AliAnalysisTaskSELambdac::UserExec: MC particles branch not found!\n");
 508       return;
 509     }
 510
 511   // load MC header
 512     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 513     if(!mcHeader) {
 514     printf("AliAnalysisTaskSELambdac::UserExec: MC header branch not found!\n");
 515     return;
 516     }
 517   }
 518
 519   Int_t n3Prong = array3Prong->GetEntriesFast();
 520
 521
 522   Int_t nOS=0;
 523   Int_t index;
 524   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 525     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 526
 527
 528     Bool_t unsetvtx=kFALSE;
 529     if(!d->GetOwnPrimaryVtx()){
 530       d->SetOwnPrimaryVtx(vtx1);
 531       unsetvtx=kTRUE;
 532     }
 533
 534     if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate)) {
 535       Int_t iPtBin = -1;
 536       Double_t ptCand = d->Pt();
 537
 538       for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
 539         if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
 540       }
 541
 542       Bool_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate);
 543       Int_t labDp=-1;
 544       Float_t deltaPx=0.;
 545       Float_t deltaPy=0.;
 546       Float_t deltaPz=0.;
 547       Float_t truePt=0.;
 548       Float_t xDecay=0.;
 549       Float_t yDecay=0.;
 550       Float_t zDecay=0.;
 551       Float_t pdgCode=-2;
 552       if(fReadMC){
 553         labDp = MatchToMCLambdac(d,arrayMC);
 554         //
 555         if(labDp>0){
 556           AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
 557           AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
 558           deltaPx=partDp->Px()-d->Px();
 559           deltaPy=partDp->Py()-d->Py();
 560           deltaPz=partDp->Pz()-d->Pz();
 561           truePt=partDp->Pt();
 562           xDecay=dg0->Xv();
 563           yDecay=dg0->Yv();
 564           zDecay=dg0->Zv();
 565           pdgCode=TMath::Abs(partDp->GetPdgCode());
 566         }else{
 567           pdgCode=-1;
 568         }
 569       }
 570
 571       Double_t invMasspKpi=-1.;
 572       Double_t invMasspiKp=-1.;
 573       Int_t pdgs[3]={0,0,0};
 574       Double_t field=aod->GetMagneticField();
 575       //apply MC PID
 576       if(fReadMC && fMCPid){
 577
 578        if(IspKpiMC(d,arrayMC,pdgs)) {
 579         invMasspKpi=d->InvMassLcpKpi();
 580         if(fUseKF){
 581          pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
 582          if(!VertexingKF(d,pdgs,field)) invMasspKpi=-1.;
 583         }
 584        }
 585        if(IspiKpMC(d,arrayMC)) {
 586         invMasspiKp=d->InvMassLcpiKp();
 587         if(fUseKF){
 588          pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
 589          if(!VertexingKF(d,pdgs,field)) invMasspiKp=-1.;
 590         }
 591        }
 592       }
 593       // aplly realistic PID
 594       if(fRealPid){
 595        if(IspKpiReal(d)) {
 596         invMasspKpi=d->InvMassLcpKpi();
 597         if(fUseKF){
 598          pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
 599          if(!VertexingKF(d,pdgs,field)) invMasspKpi=-1.;
 600         }
 601        }
 602        if(IspiKpReal(d)) {
 603         invMasspiKp=d->InvMassLcpiKp();
 604         if(fUseKF){
 605          pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
 606          if(!VertexingKF(d,pdgs,field)) invMasspiKp=-1.;
 607         }
 608        }
 609       }
 610       //apply PID using resonances
 611       if(fResPid){
 612        if(IspKpiResonant(d,field)) {
 613         invMasspKpi=d->InvMassLcpKpi();
 614         if(fUseKF){
 615          pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
 616          if(!VertexingKF(d,pdgs,field)) invMasspKpi=-1.;
 617         }
 618        }
 619        if(IspiKpResonant(d,field)) {
 620         invMasspiKp=d->InvMassLcpiKp();
 621         if(fUseKF){
 622          pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
 623          if(!VertexingKF(d,pdgs,field)) invMasspiKp=-1.;
 624         }
 625        }
 626       }
 627       // no PID
 628       if(!fResPid && !fRealPid && !fMCPid){
 629        invMasspiKp=d->InvMassLcpiKp();
 630        if(fUseKF){
 631          pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
 632          if(!VertexingKF(d,pdgs,field)) invMasspiKp=-1.;
 633         }
 634        invMasspKpi=d->InvMassLcpKpi();
 635        if(fUseKF){
 636         pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
 637         if(!VertexingKF(d,pdgs,field)) invMasspKpi=-1.;
 638        }
 639       }
 640
 641       if(invMasspiKp<0. && invMasspKpi<0.) continue;
 642
 643
 644       Float_t tmp[24];
 645       if(fFillNtuple){
 646         tmp[0]=pdgCode;
 647         tmp[1]=deltaPx;
 648         tmp[2]=deltaPy;
 649         tmp[3]=deltaPz;
 650         tmp[4]=truePt;
 651         tmp[5]=xDecay;
 652         tmp[6]=yDecay;
 653         tmp[7]=zDecay;
 654         tmp[8]=d->PtProng(0);
 655         tmp[9]=d->PtProng(1);
 656         tmp[10]=d->PtProng(2);
 657         tmp[11]=d->Pt();
 658         tmp[12]=d->CosPointingAngle();
 659         tmp[13]=d->DecayLength();
 660         tmp[14]=d->Xv();
 661         tmp[15]=d->Yv();
 662         tmp[16]=d->Zv();
 663         if(invMasspiKp>0.) tmp[17]=invMasspiKp;
 664         if(invMasspKpi>0.) tmp[17]=invMasspKpi;
 665         tmp[18]=d->GetSigmaVert();
 666         tmp[19]=d->Getd0Prong(0);
 667         tmp[20]=d->Getd0Prong(1);
 668         tmp[21]=d->Getd0Prong(2);
 669         tmp[22]=d->GetDCA();
 670         tmp[23]=d->Prodd0d0();
 671         fNtupleLambdac->Fill(tmp);
 672         PostData(3,fNtupleLambdac);
 673       }
 674       Double_t dlen=d->DecayLength();
 675       Double_t cosp=d->CosPointingAngle();
 676       Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 677       Double_t dca=d->GetDCA();
 678 Double_t ptmax=0;
 679       for(Int_t i=0;i<3;i++){
 680         if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
 681       }
 682       if(iPtBin>=0){
 683
 684         index=GetHistoIndex(iPtBin);
 685         if(invMasspiKp>0. && invMasspKpi>0.){
 686         if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp,0.5);
 687         if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi,0.5);
 688         }else{
 689          if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
 690          if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
 691         }
 692         fCosPHist[index]->Fill(cosp);
 693         fDLenHist[index]->Fill(dlen);
 694         fSumd02Hist[index]->Fill(sumD02);
 695         fPtMaxHist[index]->Fill(ptmax);
 696         fDCAHist[index]->Fill(dca);
 697
 698         if(passTightCuts){
 699          if(invMasspiKp>0. && invMasspKpi>0.){
 700           if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp,0.5);
 701           if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi,0.5);
 702          }else{
 703           if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp);
 704           if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi);
 705          }
 706         }
 707
 708         if(fReadMC){
 709           if(labDp>0) {
 710             index=GetSignalHistoIndex(iPtBin);
 711             if(invMasspiKp>0. && invMasspKpi>0.){
 712              if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp,0.5);
 713              if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi,0.5);
 714             }else{
 715              if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
 716              if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
 717             }
 718             fCosPHist[index]->Fill(cosp);
 719             fDLenHist[index]->Fill(dlen);
 720             fSumd02Hist[index]->Fill(sumD02);
 721             fPtMaxHist[index]->Fill(ptmax);
 722             fDCAHist[index]->Fill(dca);
 723             if(passTightCuts){
 724              if(invMasspiKp>0. && invMasspKpi>0.){
 725               if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp,0.5);
 726               if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi,0.5);
 727              }else{
 728               if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp);
 729               if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi);
 730              }
 731             }
 732
 733           }else{
 734             index=GetBackgroundHistoIndex(iPtBin);
 735             if(invMasspiKp>0. && invMasspKpi>0.){
 736              fMassHist[index]->Fill(invMasspiKp,0.5);
 737              fMassHist[index]->Fill(invMasspKpi,0.5);
 738             }else{
 739              if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
 740              if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
 741            }
 742             fCosPHist[index]->Fill(cosp);
 743             fDLenHist[index]->Fill(dlen);
 744             fSumd02Hist[index]->Fill(sumD02);
 745             fPtMaxHist[index]->Fill(ptmax);
 746             fDCAHist[index]->Fill(dca);
 747             if(passTightCuts){
 748             if(invMasspiKp>0. && invMasspKpi>0.){
 749                fMassHistTC[index]->Fill(invMasspiKp,0.5);
 750                fMassHistTC[index]->Fill(invMasspKpi,0.5);
 751               }else{
 752                if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp);
 753                if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi);
 754               }
 755             }
 756           }
 757         }
 758       }
 759       /*
 760       //start OS analysis
 761       if(labDp<0)fHistOSbkg->Fill(d->InvMassDplus());
 762       fHistOS->Fill(d->InvMassDplus());
 763       */
 764       nOS++;
 765     }
 766     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 767   }
 768
 769   PostData(1,fOutput);
 770   return;
 771 }
 772
 773
 774
 775 //________________________________________________________________________
 776 void AliAnalysisTaskSELambdac::Terminate(Option_t */*option*/)
 777 {
 778   // Terminate analysis
 779   //
 780   if(fDebug > 1) printf("AnalysisTaskSELambdac: Terminate() \n");
 781
 782   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 783   if (!fOutput) {
 784     printf("ERROR: fOutput not available\n");
 785     return;
 786   }
 787  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
 788
 789  TString hisname;
 790  Int_t index=0;
 791
 792
 793  for(Int_t i=0;i<fNPtBins;i++){
 794     index=GetHistoIndex(i);
 795     hisname.Form("hMassPt%d",i);
 796     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 797      hisname.Form("hCosPAllPt%d",i);
 798     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 799      hisname.Form("hDLenAllPt%d",i);
 800     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 801      hisname.Form("hSumd02AllPt%d",i);
 802      fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 803      hisname.Form("hSigVertAllPt%d",i);
 804      fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 805      hisname.Form("hPtMaxAllPt%d",i);
 806      fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 807      hisname.Form("hDCAAllPt%d",i);
 808      fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 809      hisname.Form("hMassPt%dTC",i);
 810     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 811
 812     index=GetSignalHistoIndex(i);
 813     hisname.Form("hSigPt%d",i);
 814     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 815     hisname.Form("hCosPSigPt%d",i);
 816     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 817     hisname.Form("hDLenSigPt%d",i);
 818     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 819     hisname.Form("hSumd02SigPt%d",i);
 820     fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 821     hisname.Form("hSigVertSigPt%d",i);
 822     fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 823     hisname.Form("hPtMaxSigPt%d",i);
 824     fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 825     hisname.Form("hDCASigPt%d",i);
 826     fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 827
 828     hisname.Form("hSigPt%dTC",i);
 829     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 830
 831     index=GetBackgroundHistoIndex(i);
 832     hisname.Form("hBkgPt%d",i);
 833     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 834     hisname.Form("hCosPBkgPt%d",i);
 835     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 836     hisname.Form("hDLenBkgPt%d",i);
 837     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 838     hisname.Form("hSumd02BkgPt%d",i);
 839     fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 840     hisname.Form("hSigVertBkgPt%d",i);
 841     fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 842     hisname.Form("hPtMaxBkgPt%d",i);
 843     fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 844     hisname.Form("hDCABkgPt%d",i);
 845     fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 846     hisname.Form("hBkgPt%dTC",i);
 847     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 848
 849  }
 850
 851   if(fFillNtuple){
 852     fNtupleLambdac = dynamic_cast<TNtuple*>(GetOutputData(3));
 853   }
 854
 855
 856  return;
 857 }
 858
 859 //________________________________________________________________________
 860 Int_t AliAnalysisTaskSELambdac::MatchToMCLambdac(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC) const{
 861   // check if the candidate is a Lambdac decaying in pKpi or in the resonant channels
 862   Int_t lambdacLab[3]={0,0,0};
 863   Int_t pdgs[3]={0,0,0};
 864   for(Int_t i=0;i<3;i++){
 865    AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
 866    Int_t lab=daugh->GetLabel();
 867    if(lab<0) return 0;
 868    AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab);
 869    if(!part) continue;
 870    pdgs[i]=part->GetPdgCode();
 871    Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
 872    if(partPdgcode==211 || partPdgcode==321 || partPdgcode==2212){
 873         Int_t motherLabel=part->GetMother();
 874              if(motherLabel<0) return 0;
 875    AliAODMCParticle *motherPart = (AliAODMCParticle*)arrayMC->At(motherLabel);
 876    if(!motherPart) continue;
 877         Int_t motherPdg = TMath::Abs(motherPart->GetPdgCode());
 878    if(motherPdg==4122) {
 879         if(GetLambdacDaugh(motherPart,arrayMC)){lambdacLab[i]=motherLabel;continue;}
 880               }
 881    if(motherPdg==313 || motherPdg==2224 || motherPdg==3124){
 882          Int_t GmotherLabel=motherPart->GetMother();
 883                if(GmotherLabel<0) return 0;
 884          AliAODMCParticle *GmotherPart = (AliAODMCParticle*)arrayMC->At(GmotherLabel);
 885          if(!GmotherPart) continue;
 886                Int_t GmotherPdg = TMath::Abs(GmotherPart->GetPdgCode());
 887          if(GmotherPdg==4122) {
 888                 if(GetLambdacDaugh(GmotherPart,arrayMC)) {lambdacLab[i]=GmotherLabel;continue;}
 889               }
 890         }
 891      }
 892   }
 893
 894  if(lambdacLab[0]==lambdacLab[1] && lambdacLab[1]==lambdacLab[2]) {return lambdacLab[0];}
 895  return 0;
 896
 897 }
 898 //------------------------
 899 Bool_t AliAnalysisTaskSELambdac::GetLambdacDaugh(AliAODMCParticle *part,TClonesArray *arrayMC) const{
 900  // check if the particle is a lambdac and if its decay mode is the correct one
 901  Int_t numberOfLambdac=0;
 902  if(TMath::Abs(part->GetPdgCode())!=4122) return kFALSE;
 903  Int_t daugh_tmp[2];
 904  daugh_tmp[0]=part->GetDaughter(0);
 905  daugh_tmp[1]=part->GetDaughter(1);
 906  Int_t nDaugh = (Int_t)part->GetNDaughters();
 907  if(nDaugh<2) return kFALSE;
 908  if(nDaugh>3) return kFALSE;
 909  AliAODMCParticle* pdaugh1 = (AliAODMCParticle*)arrayMC->At(part->GetDaughter(0));
 910    if(!pdaugh1) {return kFALSE;}
 911  Int_t number1 = TMath::Abs(pdaugh1->GetPdgCode());
 912  AliAODMCParticle* pdaugh2 = (AliAODMCParticle*)arrayMC->At(part->GetDaughter(1));
 913    if(!pdaugh2) {return kFALSE;}
 914   Int_t number2 = TMath::Abs(pdaugh2->GetPdgCode());
 915
 916   if(nDaugh==3){
 917    Int_t thirdDaugh=part->GetDaughter(1)-1;
 918             AliAODMCParticle* pdaugh3 = (AliAODMCParticle*)arrayMC->At(thirdDaugh);
 919   Int_t number3 = TMath::Abs(pdaugh3->GetPdgCode());
 920    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)) {
 921    numberOfLambdac++;
 922    }
 923   }
 924
 925  if(nDaugh==2){
 926
 927   //Lambda resonant
 928
 929   //Lambda -> p K*0
 930   //
 931   Int_t nfiglieK=0;
 932
 933    if((number1==2212 && number2==313)){
 934      nfiglieK=pdaugh2->GetNDaughters();
 935      if(nfiglieK!=2) return kFALSE;
 936      AliAODMCParticle* pdaughK1 = (AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
 937      AliAODMCParticle* pdaughK2 = (AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
 938      if(!pdaughK1) return kFALSE;
 939      if(!pdaughK2) return kFALSE;
 940      if((TMath::Abs(pdaughK1->GetPdgCode())==211 && TMath::Abs(pdaughK2->GetPdgCode())==321) || (TMath::Abs(pdaughK1->GetPdgCode())==321 && TMath::Abs(pdaughK2->GetPdgCode())==211)) numberOfLambdac++;
 941     }
 942
 943    if((number1==313 && number2==2212)){
 944     nfiglieK=pdaugh1->GetNDaughters();
 945     if(nfiglieK!=2) return kFALSE;
 946     AliAODMCParticle* pdaughK1 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
 947     AliAODMCParticle* pdaughK2 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
 948     if(!pdaughK1) return kFALSE;
 949     if(!pdaughK2) return kFALSE;
 950      if((TMath::Abs(pdaughK1->GetPdgCode())==211 && TMath::Abs(pdaughK2->GetPdgCode())==321) || (TMath::Abs(pdaughK1->GetPdgCode())==321 && TMath::Abs(pdaughK2->GetPdgCode())==211)) numberOfLambdac++;
 951    }
 952
 953    //Lambda -> Delta++ k
 954    Int_t nfiglieDelta=0;
 955    if(number1==321 && number2==2224){
 956     nfiglieDelta=pdaugh2->GetNDaughters();
 957     if(nfiglieDelta!=2) return kFALSE;
 958     AliAODMCParticle *pdaughD1=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
 959     AliAODMCParticle *pdaughD2=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
 960     if(!pdaughD1) return kFALSE;
 961     if(!pdaughD2) return kFALSE;
 962     if((TMath::Abs(pdaughD1->GetPdgCode())==211 && TMath::Abs(pdaughD2->GetPdgCode())==2212) || (TMath::Abs(pdaughD1->GetPdgCode())==2212 && TMath::Abs(pdaughD2->GetPdgCode())==211)) numberOfLambdac++;
 963    }
 964    if(number1==2224 && number2==321){
 965     nfiglieDelta=pdaugh1->GetNDaughters();
 966     if(nfiglieDelta!=2) return kFALSE;
 967     AliAODMCParticle* pdaughD1 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
 968     AliAODMCParticle* pdaughD2 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
 969     if(!pdaughD1) return kFALSE;
 970     if(!pdaughD2) return kFALSE;
 971     if((TMath::Abs(pdaughD1->GetPdgCode())==211 && TMath::Abs(pdaughD2->GetPdgCode())==2212) || (TMath::Abs(pdaughD1->GetPdgCode())==2212 && TMath::Abs(pdaughD2->GetPdgCode())==211)) numberOfLambdac++;
 972    }
 973
 974
 975    //Lambdac -> Lambda(1520) pi
 976    Int_t nfiglieLa=0;
 977    if(number1==3124 && number2==211){
 978     nfiglieLa=pdaugh1->GetNDaughters();
 979     if(nfiglieLa!=2) return kFALSE;
 980     AliAODMCParticle *pdaughL1=(AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
 981     AliAODMCParticle *pdaughL2=(AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
 982     if(!pdaughL1) return kFALSE;
 983     if(!pdaughL2) return kFALSE;
 984     if((TMath::Abs(pdaughL1->GetPdgCode())==321 && TMath::Abs(pdaughL2->GetPdgCode())==2212) || (TMath::Abs(pdaughL1->GetPdgCode())==2212 && TMath::Abs(pdaughL2->GetPdgCode())==321)) numberOfLambdac++;
 985    }
 986    if(number1==211 && number2==3124){
 987     nfiglieLa=pdaugh2->GetNDaughters();
 988     if(nfiglieLa!=2) return kFALSE;
 989     AliAODMCParticle *pdaughL1=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
 990     AliAODMCParticle *pdaughL2=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
 991     if(!pdaughL1) return kFALSE;
 992     if(!pdaughL2) return kFALSE;
 993     if((TMath::Abs(pdaughL1->GetPdgCode())==321 && TMath::Abs(pdaughL2->GetPdgCode())==2212) || (TMath::Abs(pdaughL1->GetPdgCode())==2212 && TMath::Abs(pdaughL2->GetPdgCode())==321)) numberOfLambdac++;
 994
 995    }
 996  }
 997
 998  if(numberOfLambdac>0) {return kTRUE;}
 999          return kFALSE;
1000 }
1001 //-----------------------------
1002 Bool_t AliAnalysisTaskSELambdac::IspKpiMC(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC,Int_t *pdgs) const{
1003   // Apply MC PID
1004    Int_t lab[3]={0,0,0};//,pdgs[3]={0,0,0};
1005    for(Int_t i=0;i<3;i++){
1006     AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
1007     lab[i]=daugh->GetLabel();
1008     if(lab[i]<0) return kFALSE;
1009     AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab[i]);
1010     if(!part) return kFALSE;
1011     pdgs[i]=TMath::Abs(part->GetPdgCode());
1012    }
1013
1014    if(pdgs[0]==2212 && pdgs[1]==321 && pdgs[2]==211) return kTRUE;
1015
1016    return kFALSE;
1017 }
1018 //-----------------------------
1019 Bool_t AliAnalysisTaskSELambdac::IspiKpMC(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC) const{
1020
1021   // Apply MC PID
1022    Int_t lab[3]={0,0,0},pdgs[3]={0,0,0};
1023    for(Int_t i=0;i<3;i++){
1024     AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
1025     lab[i]=daugh->GetLabel();
1026     if(lab[i]<0) return kFALSE;
1027     AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab[i]);
1028     if(!part) return kFALSE;
1029     pdgs[i]=TMath::Abs(part->GetPdgCode());
1030    }
1031
1032    if(pdgs[2]==2212 && pdgs[1]==321 && pdgs[0]==211) {return kTRUE;}
1033
1034    return kFALSE;
1035 }
1036 //--------------------------------------
1037 Bool_t AliAnalysisTaskSELambdac::IspiKpReal(AliAODRecoDecayHF3Prong *d) const{
1038
1039   // Apply realistic PID
1040   AliAODPidHF* pidObj=new AliAODPidHF();
1041   Bool_t type[2]={kFALSE,kFALSE};
1042   for(Int_t i=0;i<3;i++){
1043 //   Bool_t pid[3]={kFALSE,kFALSE,kFALSE};
1044   AliAODTrack *track=(AliAODTrack*)d->GetDaughter(i);
1045    //test asymmetric TPC PID
1046      Double_t plim[2]={0.6,2.};
1047      Double_t sigmas[5]={2.,1.,2.,3.,0.};
1048      pidObj->SetPLimit(plim);
1049      pidObj->SetSigma(sigmas);
1050      if(i==1) type[i]=pidObj->MatchTPCTOF(track,2,3,kFALSE);
1051      if(i==2) type[i]=pidObj->MatchTPCTOF(track,2,4,kFALSE);
1052      //pidinfo=pidObj->MatchTPCTOF(track,3,3,kFALSE);
1053
1054   }
1055
1056   if(type[0] && type[1]) {delete pidObj;return kTRUE;}
1057   delete pidObj;
1058  return kFALSE;
1059 }
1060 //--------------------------------------
1061 Bool_t AliAnalysisTaskSELambdac::IspKpiReal(AliAODRecoDecayHF3Prong *d) const{
1062
1063   // Apply realistic PID
1064   AliAODPidHF* pidObj=new AliAODPidHF();
1065   Bool_t type[2]={kFALSE,kFALSE};
1066   for(Int_t i=0;i<2;i++){
1067    //Bool_t pid[3]={kFALSE,kFALSE,kFALSE};
1068    AliAODTrack *track=(AliAODTrack*)d->GetDaughter(i);
1069      Double_t plim[2]={0.6,2.};
1070      Double_t sigmas[5]={2.,1.,2.,3.,0.};
1071      pidObj->SetPLimit(plim);
1072      pidObj->SetSigma(sigmas);
1073      if(i==1) type[i]=pidObj->MatchTPCTOF(track,2,3,kFALSE);
1074      if(i==0) type[i]=pidObj->MatchTPCTOF(track,2,4,kFALSE);
1075   }
1076   if(type[0] && type[1]) {delete pidObj;return kTRUE;}
1077   delete pidObj;
1078  return kFALSE;
1079 }
1080 //------------------------------------
1081 Bool_t AliAnalysisTaskSELambdac::VertexingKF(AliAODRecoDecayHF3Prong *d,Int_t *pdgs,Double_t field) const{
1082  // apply vertexing KF
1083    Int_t iprongs[3]={0,1,2};
1084    Double_t mass[2]={0.,0.};
1085  //topological constr
1086    AliKFParticle *Lambdac=d->ApplyVertexingKF(iprongs,3,pdgs,kTRUE,field,mass);
1087    if(!Lambdac) return kFALSE;
1088   Double_t probTot=TMath::Prob(Lambdac->GetChi2(),Lambdac->GetNDF());
1089   if(probTot<fCutsKF[0]) return kFALSE;
1090
1091   if(d->PtProng(0)<1.5 && d->Getd0Prong(0)>0.02) return kFALSE;
1092     if(d->PtProng(1)<1.5 && d->Getd0Prong(1)>0.02) return kFALSE;
1093       if(d->PtProng(2)<1.5 && d->Getd0Prong(2)>0.02) return kFALSE;
1094  //mass constr for K*
1095    Int_t ipRes[2];
1096    Int_t pdgres[2];
1097    mass[0]=0.8961;mass[1]=0.03;
1098    if(TMath::Abs(pdgs[0])==211){
1099     ipRes[0]=0;ipRes[1]=1;
1100     pdgres[0]=pdgs[0];pdgres[1]=321;
1101    }
1102    if(TMath::Abs(pdgs[2])==211){
1103     ipRes[0]=2;ipRes[1]=1;
1104     pdgres[0]=pdgs[2];pdgres[1]=321;
1105    }
1106    AliKFParticle *KappaStar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1107
1108   Double_t probKstar=TMath::Prob(KappaStar->GetChi2(),KappaStar->GetNDF());
1109   if(probKstar>fCutsKF[1]) {
1110     AliAODTrack *esdProng1=(AliAODTrack*)d->GetDaughter(ipRes[0]);
1111     AliAODTrack *esdProng2=(AliAODTrack*)d->GetDaughter(ipRes[1]);
1112     AliKFParticle prong1(*esdProng1,pdgres[0]);
1113     AliKFParticle prong2(*esdProng2,pdgres[1]);
1114     if(KappaStar->GetPt()<fCutsKF[2] && prong1.GetAngle(prong2)>fCutsKF[3]) return kFALSE;
1115   }
1116  //mass constr for Lambda
1117    mass[0]=1.520;mass[1]=0.005;
1118    if(TMath::Abs(pdgs[0])==2212){
1119     ipRes[0]=0;ipRes[1]=1;
1120     pdgres[0]=pdgs[0];pdgres[1]=pdgs[1];
1121    }
1122    if(TMath::Abs(pdgs[2])==2212){
1123     ipRes[0]=2;ipRes[1]=1;
1124     pdgres[0]=pdgs[2];pdgres[1]=pdgs[1];
1125    }
1126    AliKFParticle *Lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1127   Double_t probLa=TMath::Prob(Lambda1520->GetChi2(),Lambda1520->GetNDF());
1128   if(probLa>fCutsKF[4]) {
1129     AliAODTrack *esdProng1=(AliAODTrack*)d->GetDaughter(ipRes[0]);
1130     AliAODTrack *esdProng2=(AliAODTrack*)d->GetDaughter(ipRes[1]);
1131     AliKFParticle prong1(*esdProng1,pdgres[0]);
1132     AliKFParticle prong2(*esdProng2,pdgres[1]);
1133     if(Lambda1520->GetPt()<fCutsKF[5] && prong1.GetAngle(prong2)>fCutsKF[6]) return kFALSE;
1134   }
1135  //mass constr for Delta
1136    mass[0]=1.2;mass[1]=0.15;
1137    ipRes[0]=0;ipRes[1]=2;
1138    pdgres[0]=pdgs[0];pdgres[2]=pdgs[2];
1139    AliKFParticle *Delta=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1140   Double_t probDelta=TMath::Prob(Delta->GetChi2(),Delta->GetNDF());
1141   if(probDelta>fCutsKF[7]) {
1142     AliAODTrack *esdProng1=(AliAODTrack*)d->GetDaughter(ipRes[0]);
1143     AliAODTrack *esdProng2=(AliAODTrack*)d->GetDaughter(ipRes[1]);
1144     AliKFParticle prong1(*esdProng1,pdgres[0]);
1145     AliKFParticle prong2(*esdProng2,pdgres[1]);
1146     if(Delta->GetPt()<fCutsKF[8] && prong1.GetAngle(prong2)>fCutsKF[9]) return kFALSE;
1147   }
1148  return kTRUE;
1149 }
1150 //-------------------------------------
1151 Bool_t AliAnalysisTaskSELambdac::IspiKpResonant(AliAODRecoDecayHF3Prong *d,Double_t field) const{
1152
1153  // apply PID using the resonant channels
1154  //if lambda* -> pk
1155         Double_t mass[2]={1.520,0.005};
1156         Int_t ipRes[2]={1,2};
1157         Int_t pdgres[2]={321,2212};
1158         AliKFParticle *Lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1159         Double_t probLa=TMath::Prob(Lambda1520->GetChi2(),Lambda1520->GetNDF());
1160         if(probLa>0.1) return kTRUE;
1161  //K* -> kpi
1162 //        mass[0]=0.8961;mass[1]=0.03;
1163 //        ipRes[0]=0;ipRes[1]=1;
1164 //        pdgres[0]=211;pdgres[1]=321;
1165 //        AliKFParticle *Kstar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1166 //      Double_t probKa=TMath::Prob(Kstar->GetChi2(),Kstar->GetNDF());
1167 //      if(probKa>0.1) return kTRUE;
1168
1169  return kFALSE;
1170
1171 }
1172 //-------------------------------------
1173 Bool_t AliAnalysisTaskSELambdac::IspKpiResonant(AliAODRecoDecayHF3Prong *d,Double_t field) const{
1174
1175  // apply PID using the resonant channels
1176  //if lambda* -> pk
1177         Double_t mass[2]={1.520,0.005};
1178         Int_t ipRes[2]={0,1};
1179         Int_t pdgres[2]={2212,321};
1180         AliKFParticle *Lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1181         Double_t probLa=TMath::Prob(Lambda1520->GetChi2(),Lambda1520->GetNDF());
1182         if(probLa>0.1) return kTRUE;
1183  //K* -> kpi
1184 //        mass[0]=0.8961;mass[1]=0.03;
1185 //        ipRes[0]=1;ipRes[1]=2;
1186 //        pdgres[1]=211;pdgres[0]=321;
1187 //        AliKFParticle *Kstar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1188 //      Double_t probKa=TMath::Prob(Kstar->GetChi2(),Kstar->GetNDF());
1189 //      if(probKa>0.1) return kTRUE;
1190
1191  return kFALSE;
1192
1193 }