PWGHF/vertexingHF/AliAnalysisTaskSELambdac.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliAnalysisTaskSELambdac.cxx 63644 2013-07-23 09:04:11Z zconesa $ */
  17
  18 /////////////////////////////////////////////////////////////
  19 //
  20 // AliAnalysisTaskSE for the extraction of signal(e.g Lambdac) of heavy flavor
  21 // decay candidates with the MC truth.
  22 // Authors: r.romita@gsi.de
  23 /////////////////////////////////////////////////////////////
  24
  25 #include <TClonesArray.h>
  26 #include <TNtuple.h>
  27 #include <TCanvas.h>
  28 #include <TList.h>
  29 #include <TString.h>
  30 #include <TH1F.h>
  31 #include <TH2F.h>
  32 #include <TDatabasePDG.h>
  33
  34 #include <AliAnalysisDataSlot.h>
  35 #include <AliAnalysisDataContainer.h>
  36 #include "AliAnalysisManager.h"
  37 #include "AliAODHandler.h"
  38 #include "AliAODEvent.h"
  39 #include "AliAODVertex.h"
  40 #include "AliAODTrack.h"
  41 #include "AliAODMCHeader.h"
  42 #include "AliAODMCParticle.h"
  43 #include "AliAODRecoDecayHF3Prong.h"
  44 #include "AliAnalysisVertexingHF.h"
  45 #include "AliAnalysisTaskSE.h"
  46 #include "AliAnalysisTaskSELambdac.h"
  47 #include "AliKFParticle.h"
  48 #include "AliAODPidHF.h"
  49 #include "AliRDHFCutsLctopKpi.h"
  50 #include "AliRDHFCuts.h"
  51 #include "AliKFVertex.h"
  52 #include "AliESDVertex.h"
  53 //#include "AliAODpidUtil.h"
  54 #include "AliAODPid.h"
  55 #include "AliInputEventHandler.h"
  56 #include "AliPID.h"
  57 #include "AliNormalizationCounter.h"
  58
  59 ClassImp(AliAnalysisTaskSELambdac)
  60
  61
  62 //________________________________________________________________________
  63   AliAnalysisTaskSELambdac::AliAnalysisTaskSELambdac():
  64     AliAnalysisTaskSE(),
  65     fOutput(0),
  66     fHistNEvents(0),
  67     fhChi2(0),
  68     fhMassPtGreater3(0),
  69     fhMassPtGreater3TC(0),
  70     fhMassPtGreater3Kp(0),
  71     fhMassPtGreater3KpTC(0),
  72     fhMassPtGreater3Lpi(0),
  73     fhMassPtGreater3LpiTC(0),
  74     fhMassPtGreater3Dk(0),
  75     fhMassPtGreater3DkTC(0),
  76     fhMassPtGreater33Pr(0),
  77     fhMassPtGreater33PrTC(0),
  78     fhMassPtGreater2(0),
  79     fhMassPtGreater2TC(0),
  80     fhMassPtGreater2Kp(0),
  81     fhMassPtGreater2KpTC(0),
  82     fhMassPtGreater2Lpi(0),
  83     fhMassPtGreater2LpiTC(0),
  84     fhMassPtGreater2Dk(0),
  85     fhMassPtGreater2DkTC(0),
  86     fhMassPtGreater23Pr(0),
  87     fhMassPtGreater23PrTC(0),
  88     fNtupleLambdac(0),
  89     fUpmasslimit(2.486),
  90     fLowmasslimit(2.086),
  91     fNPtBins(0),
  92     fRDCutsAnalysis(0),
  93     fRDCutsProduction(0),
  94     fListCuts(0),
  95     fFillNtuple(kFALSE),
  96     fReadMC(kFALSE),
  97     fMCPid(kFALSE),
  98     fRealPid(kFALSE),
  99     fResPid(kTRUE),
 100     fUseKF(kFALSE),
 101     fAnalysis(kFALSE),
 102     fVHF(0),
 103     fFillVarHists(kFALSE),
 104     fMultiplicityHists(kFALSE),
 105     fPriorsHists(kFALSE),
 106     fNentries(0),
 107     fOutputMC(0),
 108     fAPriori(0),
 109     fMultiplicity(0),
 110     //fUtilPid(0),
 111     fPIDResponse(0),
 112     fCounter(0)
 113 {
 114   // Default constructor
 115   Float_t ptlims[7]={0.,2.,4.,6.,8.,12.,24.};
 116   SetPtBinLimit(7,ptlims);
 117    for(Int_t icut=0; icut<2; icut++) fCutsKF[icut]=0.;
 118    for(Int_t j=0; j<3*kMaxPtBins; j++){
 119     fMassHist[j]=0x0;
 120     fMassHistTC[j]=0x0;
 121     fMassHistLpi[j]=0x0;
 122     fMassHistLpiTC[j]=0x0;
 123     fMassHistKp[j]=0x0;
 124     fMassHistKpTC[j]=0x0;
 125     fMassHistDk[j]=0x0;
 126     fMassHistDkTC[j]=0x0;
 127     fMassHist3Pr[j]=0x0;
 128     fMassHist3PrTC[j]=0x0;
 129    }
 130 }
 131
 132 //________________________________________________________________________
 133 AliAnalysisTaskSELambdac::AliAnalysisTaskSELambdac(const char *name,Bool_t fillNtuple,AliRDHFCutsLctopKpi *lccutsana,AliRDHFCutsLctopKpi *lccutsprod):
 134   AliAnalysisTaskSE(name),
 135   fOutput(0),
 136   fHistNEvents(0),
 137   fhChi2(0),
 138   fhMassPtGreater3(0),
 139   fhMassPtGreater3TC(0),
 140   fhMassPtGreater3Kp(0),
 141   fhMassPtGreater3KpTC(0),
 142   fhMassPtGreater3Lpi(0),
 143   fhMassPtGreater3LpiTC(0),
 144   fhMassPtGreater3Dk(0),
 145   fhMassPtGreater3DkTC(0),
 146   fhMassPtGreater33Pr(0),
 147   fhMassPtGreater33PrTC(0),
 148   fhMassPtGreater2(0),
 149   fhMassPtGreater2TC(0),
 150   fhMassPtGreater2Kp(0),
 151   fhMassPtGreater2KpTC(0),
 152   fhMassPtGreater2Lpi(0),
 153   fhMassPtGreater2LpiTC(0),
 154   fhMassPtGreater2Dk(0),
 155   fhMassPtGreater2DkTC(0),
 156   fhMassPtGreater23Pr(0),
 157   fhMassPtGreater23PrTC(0),
 158   fNtupleLambdac(0),
 159   fUpmasslimit(2.486),
 160   fLowmasslimit(2.086),
 161   fNPtBins(0),
 162   fRDCutsAnalysis(lccutsana),
 163   fRDCutsProduction(lccutsprod),
 164   fListCuts(0),
 165   fFillNtuple(fillNtuple),
 166   fReadMC(kFALSE),
 167   fMCPid(kFALSE),
 168   fRealPid(kTRUE),
 169   fResPid(kFALSE),
 170   fUseKF(kFALSE),
 171   fAnalysis(kFALSE),
 172   fVHF(0),
 173   fFillVarHists(kFALSE),
 174   fMultiplicityHists(kFALSE),
 175   fPriorsHists(kFALSE),
 176   fNentries(0),
 177   fOutputMC(0),
 178   fAPriori(0),
 179   fMultiplicity(0),
 180   //fUtilPid(0),
 181   fPIDResponse(0),
 182   fCounter(0)
 183 {
 184   SetPtBinLimit(fRDCutsAnalysis->GetNPtBins()+1,
 185                 fRDCutsAnalysis->GetPtBinLimits());
 186    for(Int_t icut=0; icut<2; icut++) fCutsKF[icut]=0.;
 187    for(Int_t j=0; j<3*kMaxPtBins; j++){
 188     fMassHist[j]=0x0;
 189     fMassHistTC[j]=0x0;
 190     fMassHistLpi[j]=0x0;
 191     fMassHistLpiTC[j]=0x0;
 192     fMassHistKp[j]=0x0;
 193     fMassHistKpTC[j]=0x0;
 194     fMassHistDk[j]=0x0;
 195     fMassHistDkTC[j]=0x0;
 196     fMassHist3Pr[j]=0x0;
 197     fMassHist3PrTC[j]=0x0;
 198    }
 199   // Default constructor
 200   // Output slot #1 writes into a TList container
 201   DefineOutput(1,TList::Class());  //My private output
 202   DefineOutput(2,TList::Class());
 203   DefineOutput(3,TList::Class());
 204   DefineOutput(4,TH1F::Class());
 205   DefineOutput(5,TList::Class());
 206   DefineOutput(6,TList::Class());
 207   DefineOutput(7,AliNormalizationCounter::Class());
 208   if (fFillNtuple) {
 209     // Output slot #2 writes into a TNtuple container
 210     DefineOutput(8,TNtuple::Class());  //My private output
 211   }
 212 }
 213
 214 //________________________________________________________________________
 215 AliAnalysisTaskSELambdac::~AliAnalysisTaskSELambdac()
 216 {
 217   // Destructor
 218   if (fOutput) {
 219     delete fOutput;
 220     fOutput = 0;
 221   }
 222   if (fOutputMC) {
 223     delete fOutputMC;
 224     fOutputMC = 0;
 225   }
 226   if (fAPriori) {
 227     delete fAPriori;
 228     fAPriori = 0;
 229   }
 230   if (fMultiplicity) {
 231     delete fMultiplicity;
 232     fMultiplicity = 0;
 233   }
 234
 235   if (fVHF) {
 236     delete fVHF;
 237     fVHF = 0;
 238   }
 239
 240   if(fRDCutsAnalysis){
 241     delete fRDCutsAnalysis;
 242     fRDCutsAnalysis = 0;
 243   }
 244   if(fRDCutsProduction){
 245     delete fRDCutsProduction;
 246     fRDCutsProduction = 0;
 247   }
 248
 249   if (fListCuts) {
 250     delete fListCuts;
 251     fListCuts = 0;
 252   }
 253   if (fNentries){
 254     delete fNentries;
 255     fNentries = 0;
 256   }
 257   /*
 258     if (fUtilPid){
 259     delete fUtilPid;
 260     fUtilPid = 0;
 261     }
 262   */
 263   if (fPIDResponse) {
 264     delete  fPIDResponse;
 265   }
 266   if(fCounter){
 267    delete fCounter;
 268    fCounter = 0;
 269   }
 270
 271 }
 272 //_________________________________________________________________
 273 void  AliAnalysisTaskSELambdac::SetMassLimits(Float_t range){
 274   fUpmasslimit = 2.286+range;
 275   fLowmasslimit = 2.286-range;
 276 }
 277 //_________________________________________________________________
 278 void  AliAnalysisTaskSELambdac::SetMassLimits(Float_t lowlimit, Float_t uplimit){
 279   if(uplimit>lowlimit)
 280     {
 281       fUpmasslimit = lowlimit;
 282       fLowmasslimit = uplimit;
 283     }
 284 }
 285
 286
 287 //________________________________________________________________________
 288 void AliAnalysisTaskSELambdac::SetPtBinLimit(Int_t n, Float_t* lim){
 289   // define pt bins for analysis
 290   if(n>kMaxPtBins){
 291     printf("Max. number of Pt bins = %d\n",kMaxPtBins);
 292     fNPtBins=kMaxPtBins;
 293     fArrayBinLimits[0]=0.;
 294     fArrayBinLimits[1]=2.;
 295     fArrayBinLimits[2]=3.;
 296     fArrayBinLimits[3]=4.;
 297     for(Int_t i=4; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
 298   }else{
 299     fNPtBins=n-1;
 300     fArrayBinLimits[0]=lim[0];
 301     for(Int_t i=1; i<fNPtBins+1; i++)
 302       if(lim[i]>fArrayBinLimits[i-1]){
 303         fArrayBinLimits[i]=lim[i];
 304       }
 305       else {
 306         fArrayBinLimits[i]=fArrayBinLimits[i-1];
 307       }
 308     for(Int_t i=fNPtBins; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
 309   }
 310   if(fDebug > 1){
 311     printf("Number of Pt bins = %d\n",fNPtBins);
 312     for(Int_t i=0; i<fNPtBins; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fArrayBinLimits[i],fArrayBinLimits[i+1]);
 313   }
 314 }
 315 //_________________________________________________________________
 316 Double_t  AliAnalysisTaskSELambdac::GetPtBinLimit(Int_t ibin) const{
 317   if(ibin>fNPtBins)return -1;
 318   return fArrayBinLimits[ibin];
 319 }
 320
 321 //_________________________________________________________________
 322 void AliAnalysisTaskSELambdac::Init()
 323 {
 324   // Initialization
 325
 326   if (fDebug > 1) printf("AnalysisTaskSELambdac::Init() \n");
 327
 328   fListCuts=new TList();
 329   fListCuts->SetOwner();
 330
 331   fListCuts->Add(new AliRDHFCutsLctopKpi(*fRDCutsAnalysis));
 332   fListCuts->Add(new AliRDHFCutsLctopKpi(*fRDCutsProduction));
 333   PostData(2,fListCuts);
 334   return;
 335 }
 336
 337 //________________________________________________________________________
 338 void AliAnalysisTaskSELambdac::UserCreateOutputObjects()
 339 {
 340   // Create the output container
 341   //
 342   if (fDebug > 1) printf("AnalysisTaskSELambdac::UserCreateOutputObjects() \n");
 343
 344   // Several histograms are more conveniently managed in a TList
 345   fOutput = new TList();
 346   fOutput->SetOwner();
 347   fOutput->SetName("OutputHistos");
 348
 349   TString hisname;
 350   Int_t index=0;
 351   Int_t indexLS=0;
 352   for(Int_t i=0;i<fNPtBins;i++){
 353
 354     index=GetHistoIndex(i);
 355     indexLS=GetLSHistoIndex(i);
 356
 357     hisname.Form("hMassPt%d",i);
 358     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 359     fMassHist[index]->Sumw2();
 360     hisname.Form("hMassPt%dTC",i);
 361     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 362     fMassHistTC[index]->Sumw2();
 363
 364     hisname.Form("hMassPtLpi%d",i);
 365     fMassHistLpi[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 366     fMassHistLpi[index]->Sumw2();
 367     hisname.Form("hMassPtLpi%dTC",i);
 368     fMassHistLpiTC[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 369     fMassHistLpiTC[index]->Sumw2();
 370
 371     hisname.Form("hMassPtKp%d",i);
 372     fMassHistKp[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 373     fMassHistKp[index]->Sumw2();
 374     hisname.Form("hMassPtKp%dTC",i);
 375     fMassHistKpTC[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 376     fMassHistKpTC[index]->Sumw2();
 377     hisname.Form("hMassPtDk%d",i);
 378     fMassHistDk[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 379     fMassHistDk[index]->Sumw2();
 380     hisname.Form("hMassPtDk%dTC",i);
 381     fMassHistDkTC[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 382     fMassHistDkTC[index]->Sumw2();
 383
 384     hisname.Form("hMassPt3Pr%d",i);
 385     fMassHist3Pr[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 386     fMassHist3Pr[index]->Sumw2();
 387     hisname.Form("hMassPt3Pr%dTC",i);
 388     fMassHist3PrTC[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 389     fMassHist3PrTC[index]->Sumw2();
 390    //signal
 391     index=GetSignalHistoIndex(i);
 392     hisname.Form("hSigPt%d",i);
 393     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 394     fMassHist[index]->Sumw2();
 395     hisname.Form("hSigPt%dTC",i);
 396     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 397     fMassHistTC[index]->Sumw2();
 398     hisname.Form("hSigPtLpi%d",i);
 399     fMassHistLpi[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 400     fMassHistLpi[index]->Sumw2();
 401     hisname.Form("hSigPtLpi%dTC",i);
 402     fMassHistLpiTC[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 403     fMassHistLpiTC[index]->Sumw2();
 404
 405     hisname.Form("hSigPtKp%d",i);
 406     fMassHistKp[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 407     fMassHistKp[index]->Sumw2();
 408     hisname.Form("hSigPtKp%dTC",i);
 409     fMassHistKpTC[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 410     fMassHistKpTC[index]->Sumw2();
 411
 412     hisname.Form("hSigPtDk%d",i);
 413     fMassHistDk[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 414     fMassHistDk[index]->Sumw2();
 415     hisname.Form("hSigPtDk%dTC",i);
 416     fMassHistDkTC[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 417     fMassHistDkTC[index]->Sumw2();
 418
 419     hisname.Form("hSigPt3Pr%d",i);
 420     fMassHist3Pr[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 421     fMassHist3Pr[index]->Sumw2();
 422     hisname.Form("hSigPt3Pr%dTC",i);
 423     fMassHist3PrTC[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 424     fMassHist3PrTC[index]->Sumw2();
 425
 426     index=GetBackgroundHistoIndex(i);
 427     hisname.Form("hBkgPt%d",i);
 428     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 429     fMassHist[index]->Sumw2();
 430     hisname.Form("hBkgPt%dTC",i);
 431     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 432     fMassHistTC[index]->Sumw2();
 433     hisname.Form("hBkgPtLpi%d",i);
 434     fMassHistLpi[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 435     fMassHistLpi[index]->Sumw2();
 436     hisname.Form("hBkgPtLpi%dTC",i);
 437     fMassHistLpiTC[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 438     fMassHistLpiTC[index]->Sumw2();
 439
 440     hisname.Form("hBkgPtKp%d",i);
 441     fMassHistKp[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 442     fMassHistKp[index]->Sumw2();
 443     hisname.Form("hBkgPtKp%dTC",i);
 444     fMassHistKpTC[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 445     fMassHistKpTC[index]->Sumw2();
 446
 447     hisname.Form("hBkgPtDk%d",i);
 448     fMassHistDk[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 449     fMassHistDk[index]->Sumw2();
 450     hisname.Form("hBkgPtDk%dTC",i);
 451     fMassHistDkTC[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 452     fMassHistDkTC[index]->Sumw2();
 453
 454     hisname.Form("hBkgPt3Pr%d",i);
 455     fMassHist3Pr[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 456     fMassHist3Pr[index]->Sumw2();
 457     hisname.Form("hBkgPt3Pr%dTC",i);
 458     fMassHist3PrTC[index]=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 459     fMassHist3PrTC[index]->Sumw2();
 460   }
 461
 462   for(Int_t i=0; i<3*fNPtBins; i++){
 463     fOutput->Add(fMassHist[i]);
 464     fOutput->Add(fMassHistTC[i]);
 465     fOutput->Add(fMassHistLpi[i]);
 466     fOutput->Add(fMassHistLpiTC[i]);
 467     fOutput->Add(fMassHistKp[i]);
 468     fOutput->Add(fMassHistKpTC[i]);
 469     fOutput->Add(fMassHistDk[i]);
 470     fOutput->Add(fMassHistDkTC[i]);
 471     fOutput->Add(fMassHist3Pr[i]);
 472     fOutput->Add(fMassHist3PrTC[i]);
 473   }
 474
 475   fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
 476   fHistNEvents->Sumw2();
 477   fHistNEvents->SetMinimum(0);
 478   fOutput->Add(fHistNEvents);
 479
 480   fhChi2 = new TH1F("fhChi2", "Chi2",100,0.,10.);
 481   fhChi2->Sumw2();
 482   fOutput->Add(fhChi2);
 483
 484   fhMassPtGreater3=new TH1F("fhMassPtGreater3","Pt > 3 GeV/c",200,fLowmasslimit,fUpmasslimit);
 485   fhMassPtGreater3->Sumw2();
 486   fOutput->Add(fhMassPtGreater3);
 487   fhMassPtGreater3TC=new TH1F("fhMassPtGreater3TC","Pt > 3 GeV/c",200,fLowmasslimit,fUpmasslimit);
 488   fhMassPtGreater3TC->Sumw2();
 489   fOutput->Add(fhMassPtGreater3TC);
 490   fhMassPtGreater3Kp=new TH1F("fhMassPtGreater3Kp","Pt > 3 GeV/c",200,fLowmasslimit,fUpmasslimit);
 491   fhMassPtGreater3Kp->Sumw2();
 492   fOutput->Add(fhMassPtGreater3Kp);
 493   fhMassPtGreater3KpTC=new TH1F("fhMassPtGreater3KpTC","Pt > 3 GeV/c",200,fLowmasslimit,fUpmasslimit);
 494   fhMassPtGreater3KpTC->Sumw2();
 495   fOutput->Add(fhMassPtGreater3KpTC);
 496   fhMassPtGreater3Lpi=new TH1F("fhMassPtGreater3Lpi","Pt > 3 GeV/c",200,fLowmasslimit,fUpmasslimit);
 497   fhMassPtGreater3Lpi->Sumw2();
 498   fOutput->Add(fhMassPtGreater3Lpi);
 499   fhMassPtGreater3LpiTC=new TH1F("fhMassPtGreater3LpiTC","Pt > 3 GeV/c",200,fLowmasslimit,fUpmasslimit);
 500   fhMassPtGreater3LpiTC->Sumw2();
 501   fOutput->Add(fhMassPtGreater3LpiTC);
 502   fhMassPtGreater3Dk=new TH1F("fhMassPtGreater3Dk","Pt > 3 GeV/c",200,fLowmasslimit,fUpmasslimit);
 503   fhMassPtGreater3Dk->Sumw2();
 504   fOutput->Add(fhMassPtGreater3Dk);
 505   fhMassPtGreater3DkTC=new TH1F("fhMassPtGreater3DkTC","Pt > 3 GeV/c",200,fLowmasslimit,fUpmasslimit);
 506   fhMassPtGreater3DkTC->Sumw2();
 507   fOutput->Add(fhMassPtGreater3DkTC);
 508
 509   fhMassPtGreater33Pr=new TH1F("fhMassPtGreater33Pr","Pt > 3 GeV/c",200,fLowmasslimit,fUpmasslimit);
 510   fhMassPtGreater33Pr->Sumw2();
 511   fOutput->Add(fhMassPtGreater33Pr);
 512   fhMassPtGreater33PrTC=new TH1F("fhMassPtGreater33PrTC","Pt > 3 GeV/c",200,fLowmasslimit,fUpmasslimit);
 513   fhMassPtGreater33PrTC->Sumw2();
 514   fOutput->Add(fhMassPtGreater33PrTC);
 515
 516
 517   fhMassPtGreater2=new TH1F("fhMassPtGreater2","Pt > 2 GeV/c",200,fLowmasslimit,fUpmasslimit);
 518   fhMassPtGreater2->Sumw2();
 519   fOutput->Add(fhMassPtGreater2);
 520   fhMassPtGreater2TC=new TH1F("fhMassPtGreater2TC","Pt > 2 GeV/c",200,fLowmasslimit,fUpmasslimit);
 521   fhMassPtGreater2TC->Sumw2();
 522   fOutput->Add(fhMassPtGreater2TC);
 523   fhMassPtGreater2Kp=new TH1F("fhMassPtGreater2Kp","Pt > 2 GeV/c",200,fLowmasslimit,fUpmasslimit);
 524   fhMassPtGreater2Kp->Sumw2();
 525   fOutput->Add(fhMassPtGreater2Kp);
 526   fhMassPtGreater2KpTC=new TH1F("fhMassPtGreater2KpTC","Pt > 2 GeV/c",200,fLowmasslimit,fUpmasslimit);
 527   fhMassPtGreater2KpTC->Sumw2();
 528   fOutput->Add(fhMassPtGreater2KpTC);
 529   fhMassPtGreater2Lpi=new TH1F("fhMassPtGreater2Lpi","Pt > 2 GeV/c",200,fLowmasslimit,fUpmasslimit);
 530   fhMassPtGreater2Lpi->Sumw2();
 531   fOutput->Add(fhMassPtGreater2Lpi);
 532   fhMassPtGreater2LpiTC=new TH1F("fhMassPtGreater2LpiTC","Pt > 2 GeV/c",200,fLowmasslimit,fUpmasslimit);
 533   fhMassPtGreater2LpiTC->Sumw2();
 534   fOutput->Add(fhMassPtGreater2LpiTC);
 535   fhMassPtGreater2Dk=new TH1F("fhMassPtGreater2Dk","Pt > 2 GeV/c",200,fLowmasslimit,fUpmasslimit);
 536   fhMassPtGreater2Dk->Sumw2();
 537   fOutput->Add(fhMassPtGreater2Dk);
 538   fhMassPtGreater2DkTC=new TH1F("fhMassPtGreater2DkTC","Pt > 2 GeV/c",200,fLowmasslimit,fUpmasslimit);
 539   fhMassPtGreater2DkTC->Sumw2();
 540   fOutput->Add(fhMassPtGreater2DkTC);
 541   fhMassPtGreater23Pr=new TH1F("fhMassPtGreater23Pr","Pt > 2 GeV/c",200,fLowmasslimit,fUpmasslimit);
 542   fhMassPtGreater23Pr->Sumw2();
 543   fOutput->Add(fhMassPtGreater23Pr);
 544   fhMassPtGreater23PrTC=new TH1F("fhMassPtGreater23PrTC","Pt > 2 GeV/c",200,fLowmasslimit,fUpmasslimit);
 545   fhMassPtGreater23PrTC->Sumw2();
 546   fOutput->Add(fhMassPtGreater23PrTC);
 547
 548   fOutputMC = new TList();
 549   fOutputMC->SetOwner();
 550   fOutputMC->SetName("QAMCHistos");
 551
 552   //  const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
 553
 554   fNentries=new TH1F("fNentries", "Integral(1,2) = number of AODs *** Integral(2,3) = number of candidates selected with cuts *** Integral(3,4) = number of Lc selected with cuts *** Integral(4,5) = events with good vertex ***  Integral(5,6) = pt out of bounds", 11,-0.5,10.5);
 555
 556   //ROS: qui il bin assignment e' modellato su D0 ma sicuramente iv arie
 557   fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 558   fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
 559   fNentries->GetXaxis()->SetBinLabel(3,"nLcSelected");
 560   fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtxS");
 561   fNentries->GetXaxis()->SetBinLabel(5,"ptbin = -1");
 562   fNentries->GetXaxis()->SetBinLabel(6,"no daughter");
 563   fNentries->GetXaxis()->SetBinLabel(7,"nCandSel(Tr)");
 564   fNentries->GetXaxis()->SetBinLabel(8,"PID=0");
 565   fNentries->GetXaxis()->SetBinLabel(9,"PID=1");
 566   fNentries->GetXaxis()->SetBinLabel(10,"PID=2");
 567   fNentries->GetXaxis()->SetBinLabel(11,"PID=3");
 568   fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
 569
 570   hisname.Form("hMass");
 571   TH1F *hMassInv=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 572   fOutputMC->Add(hMassInv);
 573   hisname.Form("hbMass");
 574   TH1F *hBMassInv=new TH1F(hisname.Data(),hisname.Data(),200,fLowmasslimit,fUpmasslimit);
 575   fOutputMC->Add(hBMassInv);
 576
 577   // proton specific
 578   hisname.Form("hpTOFSignal");
 579   TH1F *hProtonTOFSignal=new TH1F(hisname.Data(),hisname.Data(),100,12000.,50000.0);
 580   fOutputMC->Add(hProtonTOFSignal);
 581   hisname.Form("hbpTOFSignal");
 582   TH1F *hBProtonTOFSignal=new TH1F(hisname.Data(),hisname.Data(),100,12000.,50000.0);
 583   fOutputMC->Add(hBProtonTOFSignal);
 584
 585   hisname.Form("hpTPCSignal");
 586   TH1F *hProtonTPCSignal=new TH1F(hisname.Data(),hisname.Data(),150,0.,150.0);
 587   fOutputMC->Add(hProtonTPCSignal);
 588   hisname.Form("hbpTPCSignal");
 589   TH1F *hBProtonTPCSignal=new TH1F(hisname.Data(),hisname.Data(),150,0.,150.0);
 590   fOutputMC->Add(hBProtonTPCSignal);
 591
 592   hisname.Form("hpptProng");
 593   TH1F *hProtonPtProng=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.0);
 594   fOutputMC->Add(hProtonPtProng);
 595   hisname.Form("hbpptProng");
 596   TH1F *hBProtonPtProng=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.0);
 597   fOutputMC->Add(hBProtonPtProng);
 598
 599   hisname.Form("hpRealTot");
 600   TH1F *hProtonRealTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 601   fOutputMC->Add(hProtonRealTot);
 602   hisname.Form("hbpRealTot");
 603   TH1F *hBProtonRealTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 604   fOutputMC->Add(hBProtonRealTot);
 605   hisname.Form("hpIDTot");
 606   TH1F *hProtonIDTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 607   fOutputMC->Add(hProtonIDTot);
 608   hisname.Form("hpIDGood");
 609   TH1F *hProtonIDGood=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 610   fOutputMC->Add(hProtonIDGood);
 611   hisname.Form("hbpIDGood");
 612   TH1F *hBProtonIDGood=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 613   fOutputMC->Add(hBProtonIDGood);
 614   hisname.Form("hbpIDTot");
 615   TH1F *hBProtonIDTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 616   fOutputMC->Add(hBProtonIDTot);
 617   hisname.Form("hnopIDp");
 618   TH1F *hnoProtonIDpTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 619   fOutputMC->Add(hnoProtonIDpTot);
 620   hisname.Form("hbnopIDp");
 621   TH1F *hBnoProtonIDpTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 622   fOutputMC->Add(hBnoProtonIDpTot);
 623
 624   hisname.Form("hpd0Prong");
 625   TH1F *hProtond0Prong=new TH1F(hisname.Data(),hisname.Data(),100,-0.1,0.1);
 626   fOutputMC->Add(hProtond0Prong);
 627   hisname.Form("hbpd0Prong");
 628   TH1F *hBProtond0Prong=new TH1F(hisname.Data(),hisname.Data(),100,-0.1,0.1);
 629   fOutputMC->Add(hBProtond0Prong);
 630   hisname.Form("hbpSignalVspTOF");
 631   TH2F *hBpSignalVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,12000.,50000.0);
 632   fOutputMC->Add(hBpSignalVspTOF);
 633   hisname.Form("hbpSignalVspTPC");
 634   TH2F *hBpSignalVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,150,0.,150.0);
 635   fOutputMC->Add(hBpSignalVspTPC);
 636   hisname.Form("hpSignalVspTOF");
 637   TH2F *hpSignalVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,12000.,50000.0);
 638   fOutputMC->Add(hpSignalVspTOF);
 639   hisname.Form("hpSignalVspTPC");
 640   TH2F *hpSignalVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,150,0.,150.0);
 641   fOutputMC->Add(hpSignalVspTPC);
 642
 643   hisname.Form("hpSigmaVspTOF");
 644   TH2F *hpSigmaVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 645   fOutputMC->Add(hpSigmaVspTOF);
 646   hisname.Form("hpSigmaVspTPC");
 647   TH2F *hpSigmaVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 648   fOutputMC->Add(hpSigmaVspTPC);
 649   hisname.Form("hbpSigmaVspTOF");
 650   TH2F *hBpSigmaVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 651   fOutputMC->Add(hBpSigmaVspTOF);
 652   hisname.Form("hbpSigmaVspTPC");
 653   TH2F *hBpSigmaVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 654   fOutputMC->Add(hBpSigmaVspTPC);
 655
 656
 657   //kaon specific
 658   hisname.Form("hKTOFSignal");
 659   TH1F *hKaonTOFSignal=new TH1F(hisname.Data(),hisname.Data(),100,12000.,50000.0);
 660   fOutputMC->Add(hKaonTOFSignal);
 661   hisname.Form("hbKTOFSignal");
 662   TH1F *hBKaonTOFSignal=new TH1F(hisname.Data(),hisname.Data(),100,12000.,50000.0);
 663   fOutputMC->Add(hBKaonTOFSignal);
 664   hisname.Form("hKTPCSignal");
 665   TH1F *hKaonTPCSignal=new TH1F(hisname.Data(),hisname.Data(),150,0.,150.0);
 666   fOutputMC->Add(hKaonTPCSignal);
 667   hisname.Form("hbKTPCSignal");
 668   TH1F *hBKaonTPCSignal=new TH1F(hisname.Data(),hisname.Data(),150,0.,150.0);
 669   fOutputMC->Add(hBKaonTPCSignal);
 670
 671   hisname.Form("hKptProng");
 672   TH1F *hKaonPtProng=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.0);
 673   fOutputMC->Add(hKaonPtProng);
 674   hisname.Form("hbKptProng");
 675   TH1F *hBKaonPtProng=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.0);
 676   fOutputMC->Add(hBKaonPtProng);
 677   hisname.Form("hKRealTot");
 678   TH1F *hKaonRealTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 679   fOutputMC->Add(hKaonRealTot);
 680   hisname.Form("hbKRealTot");
 681   TH1F *hBKaonRealTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 682   fOutputMC->Add(hBKaonRealTot);
 683   hisname.Form("hKIDGood");
 684   TH1F *hKaonIDGood=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 685   fOutputMC->Add(hKaonIDGood);
 686   hisname.Form("hKIDTot");
 687   TH1F *hKaonIDTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 688   fOutputMC->Add(hKaonIDTot);
 689   hisname.Form("hbKIDGood");
 690   TH1F *hBKaonIDGood=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 691   fOutputMC->Add(hBKaonIDGood);
 692   hisname.Form("hbKIDTot");
 693   TH1F *hBKaonIDTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 694   fOutputMC->Add(hBKaonIDTot);
 695   hisname.Form("hnokIDk");
 696   TH1F *hnoKaonIDkTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 697   fOutputMC->Add(hnoKaonIDkTot);
 698   hisname.Form("hbnokIDk");
 699   TH1F *hBnoKaonIDkTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 700   fOutputMC->Add(hBnoKaonIDkTot);
 701
 702
 703   hisname.Form("hKd0Prong");
 704   TH1F *hKaond0Prong=new TH1F(hisname.Data(),hisname.Data(),100,-0.1,0.1);
 705   fOutputMC->Add(hKaond0Prong);
 706   hisname.Form("hbKd0Prong");
 707   TH1F *hBKaond0Prong=new TH1F(hisname.Data(),hisname.Data(),100,-0.1,0.1);
 708   fOutputMC->Add(hBKaond0Prong);
 709   hisname.Form("hbKSignalVspTOF");
 710   TH2F *hbKSignalVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,12000.,50000.0);
 711   fOutputMC->Add(hbKSignalVspTOF);
 712   hisname.Form("hbKSignalVspTPC");
 713   TH2F *hbKSignalVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,150,0.,150.0);
 714   fOutputMC->Add(hbKSignalVspTPC);
 715   hisname.Form("hKSignalVspTOF");
 716   TH2F *hKSignalVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,12000.,50000.0);
 717   fOutputMC->Add(hKSignalVspTOF);
 718   hisname.Form("hKSignalVspTPC");
 719   TH2F *hKSignalVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,150,0.,150.0);
 720   fOutputMC->Add(hKSignalVspTPC);
 721   hisname.Form("hKSigmaVspTOF");
 722   TH2F *hKSigmaVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 723   fOutputMC->Add(hKSigmaVspTOF);
 724   hisname.Form("hKSigmaVspTPC");
 725   TH2F *hKSigmaVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 726   fOutputMC->Add(hKSigmaVspTPC);
 727   hisname.Form("hbKSigmaVspTOF");
 728   TH2F *hBKSigmaVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 729   fOutputMC->Add(hBKSigmaVspTOF);
 730   hisname.Form("hbKSigmaVspTPC");
 731   TH2F *hBKSigmaVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 732   fOutputMC->Add(hBKSigmaVspTPC);
 733
 734
 735   // pion specific
 736   hisname.Form("hpiTOFSignal");
 737   TH1F *hPionTOFSignal=new TH1F(hisname.Data(),hisname.Data(),100,12000.,50000.0);
 738   fOutputMC->Add(hPionTOFSignal);
 739   hisname.Form("hbpiTOFSignal");
 740   TH1F *hBPionTOFSignal=new TH1F(hisname.Data(),hisname.Data(),100,12000.,50000.0);
 741   fOutputMC->Add(hBPionTOFSignal);
 742   hisname.Form("hpiTPCSignal");
 743   TH1F *hPionTPCSignal=new TH1F(hisname.Data(),hisname.Data(),100,30.,100.0);
 744   fOutputMC->Add(hPionTPCSignal);
 745   hisname.Form("hbpiTPCSignal");
 746   TH1F *hBPionTPCSignal=new TH1F(hisname.Data(),hisname.Data(),100,30.,100.0);
 747   fOutputMC->Add(hBPionTPCSignal);
 748
 749   hisname.Form("hpiptProng");
 750   TH1F *hPionPtProng=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.0);
 751   fOutputMC->Add(hPionPtProng);
 752   hisname.Form("hbpiptProng");
 753   TH1F *hBPionPtProng=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.0);
 754   fOutputMC->Add(hBPionPtProng);
 755
 756   hisname.Form("hpiRealTot");
 757   TH1F *hPionRealTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 758   fOutputMC->Add(hPionRealTot);
 759   hisname.Form("hbpiRealTot");
 760   TH1F *hBPionRealTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 761   fOutputMC->Add(hBPionRealTot);
 762   hisname.Form("hpiIDGood");
 763   TH1F *hPionIDGood=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 764   fOutputMC->Add(hPionIDGood);
 765   hisname.Form("hpiIDTot");
 766   TH1F *hPionIDTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 767   fOutputMC->Add(hPionIDTot);
 768   hisname.Form("hbpiIDTot");
 769   TH1F *hBPionIDTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 770   fOutputMC->Add(hBPionIDTot);
 771   hisname.Form("hbpiIDGood");
 772   TH1F *hBPionIDGood=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 773   fOutputMC->Add(hBPionIDGood);
 774   hisname.Form("hnopiIDpi");
 775   TH1F *hnoPionIDpiTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 776   fOutputMC->Add(hnoPionIDpiTot);
 777   hisname.Form("hbnopiIDpi");
 778   TH1F *hBnoPionIDpiTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 779   fOutputMC->Add(hBnoPionIDpiTot);
 780
 781   hisname.Form("hpid0Prong");
 782   TH1F *hPiond0Prong=new TH1F(hisname.Data(),hisname.Data(),100,-0.1,-0.1);
 783   fOutputMC->Add(hPiond0Prong);
 784   hisname.Form("hbpid0Prong");
 785   TH1F *hBPiond0Prong=new TH1F(hisname.Data(),hisname.Data(),100,-0.1,0.1);
 786   fOutputMC->Add(hBPiond0Prong);
 787
 788   hisname.Form("hpiSignalVspTOF");
 789   TH2F *hpiSignalVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,12000.,50000.0);
 790   fOutputMC->Add(hpiSignalVspTOF);
 791   hisname.Form("hpiSignalVspTPC");
 792   TH2F *hpiSignalVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,150,0.,150.0);
 793   fOutputMC->Add(hpiSignalVspTPC);
 794   hisname.Form("hbpiSignalVspTOF");
 795   TH2F *hbpiSignalVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,12000.,50000.0);
 796   fOutputMC->Add(hbpiSignalVspTOF);
 797   hisname.Form("hbpiSignalVspTPC");
 798   TH2F *hbpiSignalVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,150,0.,150.0);
 799   fOutputMC->Add(hbpiSignalVspTPC);
 800   hisname.Form("hpiSigmaVspTOF");
 801   TH2F *hpiSigmaVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 802   fOutputMC->Add(hpiSigmaVspTOF);
 803   hisname.Form("hpiSigmaVspTPC");
 804   TH2F *hpiSigmaVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 805   fOutputMC->Add(hpiSigmaVspTPC);
 806   hisname.Form("hbpiSigmaVspTOF");
 807   TH2F *hBpiSigmaVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 808   fOutputMC->Add(hBpiSigmaVspTOF);
 809   hisname.Form("hbpiSigmaVspTPC");
 810   TH2F *hBpiSigmaVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 811   fOutputMC->Add(hBpiSigmaVspTPC);
 812
 813   // other generic
 814   hisname.Form("hLcpt");
 815   TH1F *hLcPt=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.0);
 816   fOutputMC->Add(hLcPt);
 817   hisname.Form("hbLcpt");
 818   TH1F *hBLcPt=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.0);
 819   fOutputMC->Add(hBLcPt);
 820   hisname.Form("hDist12toPrim");
 821   TH1F *hDist12Prim=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.0);
 822   fOutputMC->Add(hDist12Prim);
 823   hisname.Form("hbDist12toPrim");
 824   TH1F *hBDist12Prim=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.0);
 825   fOutputMC->Add(hBDist12Prim);
 826
 827   hisname.Form("hSigmaVert");
 828   TH1F *hSigmaVert=new TH1F(hisname.Data(),hisname.Data(),60,0.,0.06);
 829   fOutputMC->Add(hSigmaVert);
 830   hisname.Form("hbSigmaVert");
 831   TH1F *hBSigmaVert=new TH1F(hisname.Data(),hisname.Data(),60,0.,0.06);
 832   fOutputMC->Add(hBSigmaVert);
 833
 834   hisname.Form("hDCAs");
 835   TH1F *hdcas=new TH1F(hisname.Data(),hisname.Data(),200,0.,0.1);
 836   fOutputMC->Add(hdcas);
 837   hisname.Form("hbDCAs");
 838   TH1F *hBdcas=new TH1F(hisname.Data(),hisname.Data(),200,0.,0.1);
 839   fOutputMC->Add(hBdcas);
 840
 841   hisname.Form("hCosPointingAngle");
 842   TH1F *hCosPointingAngle=new TH1F(hisname.Data(),hisname.Data(),40,0.,1.);
 843   fOutputMC->Add(hCosPointingAngle);
 844   hisname.Form("hbCosPointingAngle");
 845   TH1F *hBCosPointingAngle=new TH1F(hisname.Data(),hisname.Data(),40,0.,1.);
 846   fOutputMC->Add(hBCosPointingAngle);
 847
 848   hisname.Form("hDecayLength");
 849   TH1F *hDecayLength=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 850   fOutputMC->Add(hDecayLength);
 851   hisname.Form("hbDecayLength");
 852   TH1F *hBDecayLength=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 853   fOutputMC->Add(hBDecayLength);
 854
 855   hisname.Form("hSum2");
 856   TH1F *hSum2=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 857   fOutputMC->Add(hSum2);
 858   hisname.Form("hbSum2");
 859   TH1F *hBSum2=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 860   fOutputMC->Add(hBSum2);
 861
 862   hisname.Form("hptmax");
 863   TH1F *hPtMax=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 864   fOutputMC->Add(hPtMax);
 865   hisname.Form("hbptmax");
 866   TH1F *hBPtMax=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 867   fOutputMC->Add(hBPtMax);
 868
 869
 870
 871   fAPriori = new TList(); // AdC
 872   fAPriori->SetOwner(); // AdC
 873   fAPriori->SetName("APrioriMCHistos"); // AdC
 874
 875   hisname.Form("hElIn3Prong");
 876   TH1F *hElIn3Prong=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 877   fAPriori->Add(hElIn3Prong);
 878   hisname.Form("hMuIn3Prong");
 879   TH1F *hMuIn3Prong=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 880   fAPriori->Add(hMuIn3Prong);
 881   hisname.Form("hPiIn3Prong");
 882   TH1F *hPiIn3Prong=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 883   fAPriori->Add(hPiIn3Prong);
 884   hisname.Form("hKaIn3Prong");
 885   TH1F *hKaIn3Prong=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 886   fAPriori->Add(hKaIn3Prong);
 887   hisname.Form("hPrIn3Prong");
 888   TH1F *hPrIn3Prong=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 889   fAPriori->Add(hPrIn3Prong);
 890
 891   hisname.Form("hElIn3Prong1");
 892   TH1F *hElIn3Prong1=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 893   fAPriori->Add(hElIn3Prong1);
 894   hisname.Form("hMuIn3Prong1");
 895   TH1F *hMuIn3Prong1=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 896   fAPriori->Add(hMuIn3Prong1);
 897   hisname.Form("hPiIn3Prong1");
 898   TH1F *hPiIn3Prong1=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 899   fAPriori->Add(hPiIn3Prong1);
 900   hisname.Form("hKaIn3Prong1");
 901   TH1F *hKaIn3Prong1=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 902   fAPriori->Add(hKaIn3Prong1);
 903   hisname.Form("hPrIn3Prong1");
 904   TH1F *hPrIn3Prong1=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 905   fAPriori->Add(hPrIn3Prong1);
 906
 907   hisname.Form("hElIn3Prong2");
 908   TH1F *hElIn3Prong2=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 909   fAPriori->Add(hElIn3Prong2);
 910   hisname.Form("hMuIn3Prong2");
 911   TH1F *hMuIn3Prong2=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 912   fAPriori->Add(hMuIn3Prong2);
 913   hisname.Form("hPiIn3Prong2");
 914   TH1F *hPiIn3Prong2=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 915   fAPriori->Add(hPiIn3Prong2);
 916   hisname.Form("hKaIn3Prong2");
 917   TH1F *hKaIn3Prong2=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 918   fAPriori->Add(hKaIn3Prong2);
 919   hisname.Form("hPrIn3Prong2");
 920   TH1F *hPrIn3Prong2=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 921   fAPriori->Add(hPrIn3Prong2);
 922
 923   hisname.Form("hElIn3Prong3");
 924   TH1F *hElIn3Prong3=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 925   fAPriori->Add(hElIn3Prong3);
 926   hisname.Form("hMuIn3Prong3");
 927   TH1F *hMuIn3Prong3=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 928   fAPriori->Add(hMuIn3Prong3);
 929   hisname.Form("hPiIn3Prong3");
 930   TH1F *hPiIn3Prong3=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 931   fAPriori->Add(hPiIn3Prong3);
 932   hisname.Form("hKaIn3Prong3");
 933   TH1F *hKaIn3Prong3=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 934   fAPriori->Add(hKaIn3Prong3);
 935   hisname.Form("hPrIn3Prong3");
 936   TH1F *hPrIn3Prong3=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 937   fAPriori->Add(hPrIn3Prong3);
 938
 939   hisname.Form("hElIn3Prong4");
 940   TH1F *hElIn3Prong4=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 941   fAPriori->Add(hElIn3Prong4);
 942   hisname.Form("hMuIn3Prong4");
 943   TH1F *hMuIn3Prong4=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 944   fAPriori->Add(hMuIn3Prong4);
 945   hisname.Form("hPiIn3Prong4");
 946   TH1F *hPiIn3Prong4=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 947   fAPriori->Add(hPiIn3Prong4);
 948   hisname.Form("hKaIn3Prong4");
 949   TH1F *hKaIn3Prong4=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 950   fAPriori->Add(hKaIn3Prong4);
 951   hisname.Form("hPrIn3Prong4");
 952   TH1F *hPrIn3Prong4=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 953   fAPriori->Add(hPrIn3Prong4);
 954
 955   hisname.Form("hElIn3Prong5");
 956   TH1F *hElIn3Prong5=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 957   fAPriori->Add(hElIn3Prong5);
 958   hisname.Form("hMuIn3Prong5");
 959   TH1F *hMuIn3Prong5=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 960   fAPriori->Add(hMuIn3Prong5);
 961   hisname.Form("hPiIn3Prong5");
 962   TH1F *hPiIn3Prong5=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 963   fAPriori->Add(hPiIn3Prong5);
 964   hisname.Form("hKaIn3Prong5");
 965   TH1F *hKaIn3Prong5=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 966   fAPriori->Add(hKaIn3Prong5);
 967   hisname.Form("hPrIn3Prong5");
 968   TH1F *hPrIn3Prong5=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 969   fAPriori->Add(hPrIn3Prong5);
 970
 971   hisname.Form("hElIn3Prong6");
 972   TH1F *hElIn3Prong6=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 973   fAPriori->Add(hElIn3Prong6);
 974   hisname.Form("hMuIn3Prong6");
 975   TH1F *hMuIn3Prong6=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 976   fAPriori->Add(hMuIn3Prong6);
 977   hisname.Form("hPiIn3Prong6");
 978   TH1F *hPiIn3Prong6=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 979   fAPriori->Add(hPiIn3Prong6);
 980   hisname.Form("hKaIn3Prong6");
 981   TH1F *hKaIn3Prong6=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 982   fAPriori->Add(hKaIn3Prong6);
 983   hisname.Form("hPrIn3Prong6");
 984   TH1F *hPrIn3Prong6=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 985   fAPriori->Add(hPrIn3Prong6);
 986
 987
 988   fMultiplicity = new TList(); // AdC
 989   fMultiplicity->SetOwner(); // AdC
 990   fMultiplicity->SetName("MultiplicityMCHistos"); // AdC
 991
 992   hisname.Form("hLcinEvent");
 993   TH1I*hLcinEvent=new TH1I(hisname.Data(),hisname.Data(),3,-1,2);
 994   fMultiplicity->Add(hLcinEvent);
 995
 996   hisname.Form("hPrimariesvsAOD");
 997   TH2I *hPrimariesvsAOD=new TH2I(hisname.Data(),hisname.Data(),300,0,300,300,0,300);
 998   fMultiplicity->Add(hPrimariesvsAOD);
 999
1000   hisname.Form("hMultiplicityInLcEvent");
1001   TH1I *hMultiplicityInLcEvent=new TH1I(hisname.Data(),hisname.Data(),300,0,300);
1002   fMultiplicity->Add(hMultiplicityInLcEvent);
1003   hisname.Form("h2MultiplicityInLcEvent");
1004   TH1I *h2MultiplicityInLcEvent=new TH1I(hisname.Data(),hisname.Data(),300,0,300);
1005   fMultiplicity->Add(h2MultiplicityInLcEvent);
1006   hisname.Form("hAll2MultiplicityInEvent");
1007   TH1I *hAll2MultiplicityInEvent=new TH1I(hisname.Data(),hisname.Data(),300,0,300);
1008   fMultiplicity->Add(hAll2MultiplicityInEvent);
1009   hisname.Form("hAllMultiplicityInEvent");
1010   TH1I *hAllMultiplicityInEvent=new TH1I(hisname.Data(),hisname.Data(),300,0,300);
1011   fMultiplicity->Add(hAllMultiplicityInEvent);
1012   hisname.Form("hAllMultiplicityPrimaryInEvent");
1013   TH1I *hAllMultiplicityPrimaryInEvent=new TH1I(hisname.Data(),hisname.Data(),300,0,300);
1014   fMultiplicity->Add(hAllMultiplicityPrimaryInEvent);
1015   hisname.Form("hAll2MultiplicityPrimaryInEvent");
1016   TH1I *hAll2MultiplicityPrimaryInEvent=new TH1I(hisname.Data(),hisname.Data(),300,0,300);
1017   fMultiplicity->Add(hAll2MultiplicityPrimaryInEvent);
1018   hisname.Form("hMultiplicityInEvent");
1019   TH1I *hMultiplicityInEvent=new TH1I(hisname.Data(),hisname.Data(),300,0.,300);
1020   fMultiplicity->Add(hMultiplicityInEvent);
1021   hisname.Form("hMultiplicityIn3ProngLC");
1022   TH1I *hMultiplicityIn3ProngLC=new TH1I(hisname.Data(),hisname.Data(),300,0,300);
1023   fMultiplicity->Add(hMultiplicityIn3ProngLC);
1024   hisname.Form("hMultiplicityInLCpid");
1025   TH1I *hMultiplicityInLCpid=new TH1I(hisname.Data(),hisname.Data(),300,0,300);
1026   fMultiplicity->Add(hMultiplicityInLCpid);
1027   hisname.Form("hMultiplicityInLCmc");
1028   TH1I *hMultiplicityInLCmc=new TH1I(hisname.Data(),hisname.Data(),300,0,300);
1029   fMultiplicity->Add(hMultiplicityInLCmc);
1030   hisname.Form("hMultiplicityInLCNomc");
1031   TH1I *hMultiplicityInLCNomc=new TH1I(hisname.Data(),hisname.Data(),300,0,300);
1032   fMultiplicity->Add(hMultiplicityInLCNomc);
1033   hisname.Form("hMultiplicityYesC");
1034   TH1I *hMultiplicityYesC=new TH1I(hisname.Data(),hisname.Data(),300,0,300);
1035   fMultiplicity->Add(hMultiplicityYesC);
1036   hisname.Form("hMultiplicityYesB");
1037   TH1I *hMultiplicityYesB=new TH1I(hisname.Data(),hisname.Data(),300,0,300);
1038   fMultiplicity->Add(hMultiplicityYesB);
1039   hisname.Form("hMultiplicityJPsi");
1040   TH1I *hMultiplicityJPsi=new TH1I(hisname.Data(),hisname.Data(),300,0,300);
1041   fMultiplicity->Add(hMultiplicityJPsi);
1042
1043
1044   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
1045   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1046   fPIDResponse = inputHandler->GetPIDResponse();
1047
1048   if(fRDCutsProduction->GetIsUsePID()){
1049     fRDCutsProduction->GetPidHF()->SetPidResponse(fPIDResponse);
1050     fRDCutsProduction->GetPidpion()->SetPidResponse(fPIDResponse);
1051     fRDCutsProduction->GetPidprot()->SetPidResponse(fPIDResponse);
1052     //fUtilPid=new AliAODpidUtil(fPIDResponse);
1053     fRDCutsProduction->GetPidHF()->SetOldPid(kFALSE);
1054     fRDCutsProduction->GetPidpion()->SetOldPid(kFALSE);
1055     fRDCutsProduction->GetPidprot()->SetOldPid(kFALSE);
1056   }
1057   if(fRDCutsAnalysis->GetIsUsePID()){
1058     fRDCutsAnalysis->GetPidHF()->SetPidResponse(fPIDResponse);
1059     fRDCutsAnalysis->GetPidpion()->SetPidResponse(fPIDResponse);
1060     fRDCutsAnalysis->GetPidprot()->SetPidResponse(fPIDResponse);
1061     fRDCutsAnalysis->GetPidHF()->SetOldPid(kFALSE);
1062     fRDCutsAnalysis->GetPidpion()->SetOldPid(kFALSE);
1063     fRDCutsAnalysis->GetPidprot()->SetOldPid(kFALSE);
1064   }
1065
1066   PostData(1,fOutput);
1067   if (fFillVarHists) PostData(3,fOutputMC);
1068   PostData(4,fNentries);
1069   if (fPriorsHists) PostData(5,fAPriori);
1070   if (fMultiplicityHists) PostData(6,fMultiplicity);
1071    TString normName="NormalizationCounter";
1072    AliAnalysisDataContainer *cont = GetOutputSlot(7)->GetContainer();
1073   if(cont)normName=(TString)cont->GetName();
1074   fCounter = new AliNormalizationCounter(normName.Data());
1075    fCounter->Init();
1076     PostData(7,fCounter);
1077   if (fFillNtuple) {
1078     //OpenFile(3); // 2 is the slot number of the ntuple
1079     fNtupleLambdac = new TNtuple("fNtupleLambdac","D +",
1080                                  "pdg:Px:Py:Pz:PtTrue:VxTrue:VyTrue:VzTrue:Ptpi:PtK:Ptpi2:PtRec:PointingAngle:DecLeng:VxRec:VyRec:VzRec:InvMass:sigvert:d0Pi:d0K:d0Pi2:dca:d0square");
1081     PostData(8,fNtupleLambdac);
1082   }
1083
1084   return;
1085 }
1086
1087 //________________________________________________________________________
1088 void AliAnalysisTaskSELambdac::UserExec(Option_t */*option*/)
1089 {
1090   // Execute analysis for current event:
1091   // heavy flavor candidates association to MC truth
1092
1093   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
1094   //tmp
1095   fHistNEvents->Fill(0); // count event
1096   // Post the data already here
1097
1098   TClonesArray *array3Prong = 0;
1099   TClonesArray *arrayLikeSign =0;
1100   if(!aod && AODEvent() && IsStandardAOD()) {
1101     // In case there is an AOD handler writing a standard AOD, use the AOD
1102     // event in memory rather than the input (ESD) event.
1103     aod = dynamic_cast<AliAODEvent*> (AODEvent());
1104     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
1105     // have to taken from the AOD event hold by the AliAODExtension
1106     AliAODHandler* aodHandler = (AliAODHandler*)
1107       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
1108     if(aodHandler->GetExtensions()) {
1109       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
1110       AliAODEvent *aodFromExt = ext->GetAOD();
1111       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
1112       arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
1113     }
1114   } else if(aod) {
1115     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
1116     arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
1117   }
1118
1119   if(!aod) return;
1120
1121   TClonesArray *arrayMC=0;
1122   AliAODMCHeader *mcHeader=0;
1123
1124   // load MC particles
1125   if(fReadMC){
1126
1127     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1128     if(!arrayMC) {
1129       AliError("AliAnalysisTaskSELambdac::UserExec: MC particles branch not found!\n");
1130       return;
1131     }
1132
1133
1134     // load MC header
1135     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1136     if(!mcHeader) {
1137       AliError("AliAnalysisTaskSELambdac::UserExec: MC header branch not found!\n");
1138       return;
1139     }
1140   }
1141
1142   TString fillthis="";
1143   Int_t numberOfPrimaries= NumberPrimaries(aod);
1144
1145   if (fMultiplicityHists && fReadMC) {
1146     fillthis="hPrimariesvsAOD";
1147     ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(aod->GetNTracks(),numberOfPrimaries);
1148
1149     fillthis="hAll2MultiplicityInEvent";
1150     ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(aod->GetNTracks());
1151
1152     fillthis="hAll2MultiplicityPrimaryInEvent";
1153     ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1154
1155     if (IsThereAGeneratedLc(arrayMC)) {
1156       fillthis="h2MultiplicityInLcEvent";
1157       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1158     }
1159
1160   }
1161
1162   if(!array3Prong || !aod) {
1163     AliError("AliAnalysisTaskSELambdac::UserExec: Charm3Prong branch not found!\n");
1164     return;
1165   }
1166   if(!arrayLikeSign) {
1167     AliDebug(2,"AliAnalysisTaskSELambdac::UserExec: LikeSign3Prong branch not found!\n");
1168     //  return;
1169   }
1170
1171   // fix for temporary bug in ESDfilter
1172   // the AODs with null vertex pointer didn't pass the PhysSel
1173   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
1174   if(!vtx1 || TMath::Abs(aod->GetMagneticField())<0.001) return;
1175
1176   fNentries->Fill(0);
1177   fCounter->StoreEvent(aod,fRDCutsProduction,fReadMC);
1178   TString trigclass=aod->GetFiredTriggerClasses();
1179   if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fNentries->Fill(14);
1180   Bool_t isEvSelAnCuts,isEvSelProdCuts;
1181   isEvSelProdCuts=fRDCutsProduction->IsEventSelected(aod);
1182   isEvSelAnCuts=fRDCutsAnalysis->IsEventSelected(aod);
1183   if(!isEvSelProdCuts){
1184     if(fRDCutsProduction->GetWhyRejection()==1) // rejected for pileup
1185       fNentries->Fill(13);
1186     return;
1187   }
1188
1189   Bool_t isThereA3prongWithGoodTracks = kFALSE;
1190   Bool_t isThereA3ProngLcKine = kFALSE;
1191   Bool_t isThereA3ProngLcKineANDpid = kFALSE;
1192   Bool_t isThereA3ProngLcMC = kFALSE;
1193   Bool_t isThereA3ProngCyes = kFALSE;
1194   Bool_t isThereA3ProngByes = kFALSE;
1195   Bool_t isThereA3ProngJPsi = kFALSE;
1196
1197   Int_t n3Prong = array3Prong->GetEntriesFast();
1198   Int_t nSelectedloose[1]={0};
1199   Int_t nSelectedtight[1]={0};
1200
1201   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
1202     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
1203
1204     Bool_t unsetvtx=kFALSE;
1205     if(!d->GetOwnPrimaryVtx()){
1206       d->SetOwnPrimaryVtx(vtx1);
1207       unsetvtx=kTRUE;
1208     }
1209
1210     //if(d->GetSelectionMap()) {if(!d->HasSelectionBit(AliRDHFCuts::kLcCuts)) continue;}
1211     //if(d->GetSelectionMap()) {if(!d->HasSelectionBit(AliRDHFCuts::kLcPID)) continue;}
1212
1213     Int_t isSelectedTracks = fRDCutsProduction->IsSelected(d,AliRDHFCuts::kTracks,aod);
1214     if(!isSelectedTracks) continue;
1215
1216     isThereA3prongWithGoodTracks=kTRUE;
1217
1218     if (fRDCutsProduction->IsInFiducialAcceptance(d->Pt(),d->Y(4122))) {fNentries->Fill(6);}else{continue;}
1219
1220     Int_t ptbin=fRDCutsProduction->PtBin(d->Pt());
1221     if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
1222
1223     FillMassHists(aod,d,arrayMC,fRDCutsProduction,nSelectedloose,nSelectedtight);
1224     if(fFillVarHists) FillVarHists(d,arrayMC,fRDCutsProduction,/*fOutputMC,*/aod);
1225
1226     if (fPriorsHists && fReadMC)
1227       FillAPrioriConcentrations(d, fRDCutsProduction, aod, arrayMC);
1228     if (fMultiplicityHists && fReadMC)
1229       MultiplicityStudies(d, fRDCutsProduction, aod, arrayMC,
1230                           isThereA3ProngLcKine,isThereA3ProngLcKineANDpid,isThereA3ProngLcMC,
1231                           isThereA3ProngCyes,isThereA3ProngByes,isThereA3ProngJPsi);
1232
1233     /*
1234     //start OS analysis
1235     if(labDp<0)fHistOSbkg->Fill(d->InvMassDplus());
1236     fHistOS->Fill(d->InvMassDplus());
1237     */
1238
1239     if(unsetvtx) d->UnsetOwnPrimaryVtx();
1240   }
1241   fCounter->StoreCandidates(aod,nSelectedloose[0],kTRUE);
1242   fCounter->StoreCandidates(aod,nSelectedtight[0],kFALSE);
1243
1244   if (fMultiplicityHists && fReadMC) {
1245     fillthis="hAllMultiplicityInEvent";
1246     ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(aod->GetNTracks());
1247
1248     fillthis="hAllMultiplicityPrimaryInEvent";
1249     ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1250
1251     if (IsThereAGeneratedLc(arrayMC)) {
1252       fillthis="hMultiplicityInLcEvent";
1253       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1254     }
1255
1256     if (isThereA3prongWithGoodTracks) {
1257       fillthis="hMultiplicityInEvent";
1258       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1259     }
1260     if (isThereA3ProngLcKine) {
1261       fillthis="hMultiplicityIn3ProngLC";
1262       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1263     }
1264     if (isThereA3ProngLcKineANDpid) {
1265       fillthis="hMultiplicityInLCpid";
1266       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1267     }
1268     if (isThereA3ProngLcMC) {
1269       fillthis="hMultiplicityInLCmc";
1270       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1271     }
1272     if (isThereA3ProngLcKine && !isThereA3ProngLcMC) {
1273       fillthis="hMultiplicityInLCNomc";
1274       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1275     }
1276
1277     if (isThereA3ProngCyes) {
1278       fillthis="hMultiplicityYesC";
1279       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1280     }
1281     if (isThereA3ProngByes) {
1282       fillthis="hMultiplicityYesB";
1283       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1284     }
1285     if (isThereA3ProngJPsi) {
1286       fillthis="hMultiplicityJPsi";
1287       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1288     }
1289
1290   }
1291
1292
1293   PostData(1,fOutput);
1294   if (fFillVarHists) PostData(3,fOutputMC);
1295   if (fPriorsHists) PostData(5,fAPriori);
1296   if (fMultiplicityHists) PostData(6,fMultiplicity);
1297
1298   PostData(4,fNentries);
1299   PostData(7,fCounter);
1300
1301   return;
1302 }
1303
1304
1305
1306 //________________________________________________________________________
1307 void AliAnalysisTaskSELambdac::Terminate(Option_t */*option*/)
1308 {
1309   // Terminate analysis
1310   //
1311   if (fDebug > 1) printf("AnalysisTaskSELambdac: Terminate() \n");
1312
1313   fOutput = dynamic_cast<TList*> (GetOutputData(1));
1314   if (!fOutput) {
1315     AliError("ERROR: fOutput not available\n");
1316     return;
1317   }
1318   //fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
1319
1320   if(fFillNtuple){
1321     fNtupleLambdac = dynamic_cast<TNtuple*>(GetOutputData(3));
1322   }
1323
1324
1325   return;
1326 }
1327
1328 //________________________________________________________________________
1329 Int_t AliAnalysisTaskSELambdac::MatchToMCLambdac(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC) const{
1330   // check if the candidate is a Lambdac decaying in pKpi or in the resonant channels
1331   Int_t lambdacLab[3]={0,0,0};
1332   Int_t pdgs[3]={0,0,0};
1333   for(Int_t i=0;i<3;i++){
1334     AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
1335     Int_t lab=daugh->GetLabel();
1336     if(lab<0) return 0;
1337     AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab);
1338     if(!part) continue;
1339     pdgs[i]=part->GetPdgCode();
1340     Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
1341     if(partPdgcode==211 || partPdgcode==321 || partPdgcode==2212){
1342       Int_t motherLabel=part->GetMother();
1343       if(motherLabel<0) return 0;
1344       AliAODMCParticle *motherPart = (AliAODMCParticle*)arrayMC->At(motherLabel);
1345       if(!motherPart) continue;
1346       Int_t motherPdg = TMath::Abs(motherPart->GetPdgCode());
1347       if(motherPdg==4122) {
1348         if(GetLambdacDaugh(motherPart,arrayMC)){lambdacLab[i]=motherLabel;continue;}
1349       }
1350       if(motherPdg==313 || motherPdg==2224 || motherPdg==3124){
1351         Int_t granMotherLabel=motherPart->GetMother();
1352         if(granMotherLabel<0) return 0;
1353         AliAODMCParticle *granMotherPart = (AliAODMCParticle*)arrayMC->At(granMotherLabel);
1354         if(!granMotherPart) continue;
1355         Int_t granMotherPdg  = TMath::Abs(granMotherPart->GetPdgCode());
1356         if(granMotherPdg ==4122) {
1357           if(GetLambdacDaugh(granMotherPart,arrayMC)) {lambdacLab[i]=granMotherLabel;continue;}
1358         }
1359       }
1360     }
1361   }
1362
1363   if(lambdacLab[0]==lambdacLab[1] && lambdacLab[1]==lambdacLab[2]) {return lambdacLab[0];}
1364   return 0;
1365
1366 }
1367 //------------------------
1368 Bool_t AliAnalysisTaskSELambdac::GetLambdacDaugh(AliAODMCParticle *part,TClonesArray *arrayMC) const{
1369   // check if the particle is a lambdac and if its decay mode is the correct one
1370   Int_t numberOfLambdac=0;
1371   if(TMath::Abs(part->GetPdgCode())!=4122) return kFALSE;
1372   Int_t daughTmp[2];
1373   daughTmp[0]=part->GetDaughter(0);
1374   daughTmp[1]=part->GetDaughter(1);
1375   Int_t nDaugh = (Int_t)part->GetNDaughters();
1376   if(nDaugh<2) return kFALSE;
1377   if(nDaugh>3) return kFALSE;
1378   AliAODMCParticle* pdaugh1 = (AliAODMCParticle*)arrayMC->At(part->GetDaughter(0));
1379   if(!pdaugh1) {return kFALSE;}
1380   Int_t number1 = TMath::Abs(pdaugh1->GetPdgCode());
1381   AliAODMCParticle* pdaugh2 = (AliAODMCParticle*)arrayMC->At(part->GetDaughter(1));
1382   if(!pdaugh2) {return kFALSE;}
1383   Int_t number2 = TMath::Abs(pdaugh2->GetPdgCode());
1384
1385   if(nDaugh==3){
1386     Int_t thirdDaugh=part->GetDaughter(1)-1;
1387     AliAODMCParticle* pdaugh3 = (AliAODMCParticle*)arrayMC->At(thirdDaugh);
1388     Int_t number3 = TMath::Abs(pdaugh3->GetPdgCode());
1389     if((number1==321 && number2==211 && number3==2212) ||
1390        (number1==211 && number2==321 && number3==2212) ||
1391        (number1==211 && number2==2212 && number3==321) ||
1392        (number1==321 && number2==2212 && number3==211) ||
1393        (number1==2212 && number2==321 && number3==211) ||
1394        (number1==2212 && number2==211 && number3==321)) {
1395       numberOfLambdac++;
1396     }
1397   }
1398
1399   if(nDaugh==2){
1400
1401     //Lambda resonant
1402
1403     //Lambda -> p K*0
1404     //
1405     Int_t nfiglieK=0;
1406
1407     if((number1==2212 && number2==313)){
1408       nfiglieK=pdaugh2->GetNDaughters();
1409       if(nfiglieK!=2) return kFALSE;
1410       AliAODMCParticle* pdaughK1 = (AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
1411       AliAODMCParticle* pdaughK2 = (AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
1412       if(!pdaughK1) return kFALSE;
1413       if(!pdaughK2) return kFALSE;
1414       if((TMath::Abs(pdaughK1->GetPdgCode())==211 && TMath::Abs(pdaughK2->GetPdgCode())==321) || (TMath::Abs(pdaughK1->GetPdgCode())==321 && TMath::Abs(pdaughK2->GetPdgCode())==211)) numberOfLambdac++;
1415     }
1416
1417     if((number1==313 && number2==2212)){
1418       nfiglieK=pdaugh1->GetNDaughters();
1419       if(nfiglieK!=2) return kFALSE;
1420       AliAODMCParticle* pdaughK1 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
1421       AliAODMCParticle* pdaughK2 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
1422       if(!pdaughK1) return kFALSE;
1423       if(!pdaughK2) return kFALSE;
1424       if((TMath::Abs(pdaughK1->GetPdgCode())==211 && TMath::Abs(pdaughK2->GetPdgCode())==321) || (TMath::Abs(pdaughK1->GetPdgCode())==321 && TMath::Abs(pdaughK2->GetPdgCode())==211)) numberOfLambdac++;
1425     }
1426
1427     //Lambda -> Delta++ k
1428     Int_t nfiglieDelta=0;
1429     if(number1==321 && number2==2224){
1430       nfiglieDelta=pdaugh2->GetNDaughters();
1431       if(nfiglieDelta!=2) return kFALSE;
1432       AliAODMCParticle *pdaughD1=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
1433       AliAODMCParticle *pdaughD2=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
1434       if(!pdaughD1) return kFALSE;
1435       if(!pdaughD2) return kFALSE;
1436       if((TMath::Abs(pdaughD1->GetPdgCode())==211 && TMath::Abs(pdaughD2->GetPdgCode())==2212) || (TMath::Abs(pdaughD1->GetPdgCode())==2212 && TMath::Abs(pdaughD2->GetPdgCode())==211)) numberOfLambdac++;
1437     }
1438     if(number1==2224 && number2==321){
1439       nfiglieDelta=pdaugh1->GetNDaughters();
1440       if(nfiglieDelta!=2) return kFALSE;
1441       AliAODMCParticle* pdaughD1 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
1442       AliAODMCParticle* pdaughD2 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
1443       if(!pdaughD1) return kFALSE;
1444       if(!pdaughD2) return kFALSE;
1445       if((TMath::Abs(pdaughD1->GetPdgCode())==211 && TMath::Abs(pdaughD2->GetPdgCode())==2212) || (TMath::Abs(pdaughD1->GetPdgCode())==2212 && TMath::Abs(pdaughD2->GetPdgCode())==211)) numberOfLambdac++;
1446     }
1447
1448
1449     //Lambdac -> Lambda(1520) pi
1450     Int_t nfiglieLa=0;
1451     if(number1==3124 && number2==211){
1452       nfiglieLa=pdaugh1->GetNDaughters();
1453       if(nfiglieLa!=2) return kFALSE;
1454       AliAODMCParticle *pdaughL1=(AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
1455       AliAODMCParticle *pdaughL2=(AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
1456       if(!pdaughL1) return kFALSE;
1457       if(!pdaughL2) return kFALSE;
1458       if((TMath::Abs(pdaughL1->GetPdgCode())==321 && TMath::Abs(pdaughL2->GetPdgCode())==2212) || (TMath::Abs(pdaughL1->GetPdgCode())==2212 && TMath::Abs(pdaughL2->GetPdgCode())==321)) numberOfLambdac++;
1459     }
1460     if(number1==211 && number2==3124){
1461       nfiglieLa=pdaugh2->GetNDaughters();
1462       if(nfiglieLa!=2) return kFALSE;
1463       AliAODMCParticle *pdaughL1=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
1464       AliAODMCParticle *pdaughL2=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
1465       if(!pdaughL1) return kFALSE;
1466       if(!pdaughL2) return kFALSE;
1467       if((TMath::Abs(pdaughL1->GetPdgCode())==321 && TMath::Abs(pdaughL2->GetPdgCode())==2212) || (TMath::Abs(pdaughL1->GetPdgCode())==2212 && TMath::Abs(pdaughL2->GetPdgCode())==321)) numberOfLambdac++;
1468
1469     }
1470   }
1471
1472   if(numberOfLambdac>0) {return kTRUE;}
1473   return kFALSE;
1474 }
1475 //-----------------------------
1476 Bool_t AliAnalysisTaskSELambdac::IspKpiMC(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC) const{
1477   // Apply MC PID
1478   Int_t lab[3]={0,0,0},pdgs[3]={0,0,0};
1479   for(Int_t i=0;i<3;i++){
1480     AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
1481     lab[i]=daugh->GetLabel();
1482     if(lab[i]<0) return kFALSE;
1483     AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab[i]);
1484     if(!part) return kFALSE;
1485     pdgs[i]=TMath::Abs(part->GetPdgCode());
1486   }
1487
1488   if(pdgs[0]==2212 && pdgs[1]==321 && pdgs[2]==211) return kTRUE;
1489
1490   return kFALSE;
1491 }
1492 //-----------------------------
1493 Bool_t AliAnalysisTaskSELambdac::IspiKpMC(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC) const{
1494
1495   // Apply MC PID
1496   Int_t lab[3]={0,0,0},pdgs[3]={0,0,0};
1497   for(Int_t i=0;i<3;i++){
1498     AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
1499     lab[i]=daugh->GetLabel();
1500     if(lab[i]<0) return kFALSE;
1501     AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab[i]);
1502     if(!part) return kFALSE;
1503     pdgs[i]=TMath::Abs(part->GetPdgCode());
1504   }
1505
1506   if(pdgs[2]==2212 && pdgs[1]==321 && pdgs[0]==211) {return kTRUE;}
1507
1508   return kFALSE;
1509 }
1510 //--------------------------------------
1511 Bool_t AliAnalysisTaskSELambdac::VertexingKF(AliAODRecoDecayHF3Prong *d,Int_t *pdgs,Double_t field) const{
1512   // apply vertexing KF
1513   Int_t iprongs[3]={0,1,2};
1514   Double_t mass[2]={0.,0.};
1515   Bool_t constraint=kFALSE;
1516   if(fCutsKF[0]>0.)constraint=kTRUE;
1517   //topological constr
1518   AliKFParticle *lambdac=d->ApplyVertexingKF(iprongs,3,pdgs,constraint,field,mass);
1519   if(!lambdac) return kFALSE;
1520   if(lambdac->GetChi2()/lambdac->GetNDF()>fCutsKF[1]) return kFALSE;
1521   return kTRUE;
1522 }
1523 //-------------------------------------
1524 void AliAnalysisTaskSELambdac::IspiKpResonant(AliAODRecoDecayHF3Prong *d,Double_t field,Int_t *resNumber) const{
1525
1526   // apply PID using the resonant channels
1527   //if lambda* -> pk
1528   Double_t mass[2]={1.520,0.005};
1529   Int_t ipRes[2]={1,2};
1530   Int_t pdgres[2]={321,2212};
1531   AliKFParticle *lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1532   Double_t probLa=TMath::Prob(lambda1520->GetChi2(),lambda1520->GetNDF());
1533   if(probLa>0.1) resNumber[0]=1;
1534   //K* -> kpi
1535   mass[0]=0.8961;mass[1]=0.03;
1536   ipRes[0]=0;ipRes[1]=1;
1537   pdgres[0]=211;pdgres[1]=321;
1538   AliKFParticle *kstar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1539   Double_t probKa=TMath::Prob(kstar->GetChi2(),kstar->GetNDF());
1540   if(probKa>0.1) resNumber[1]=1;
1541
1542  // Delta++
1543   mass[0]=1.232;mass[1]=0.15;
1544   ipRes[0]=0;ipRes[1]=2;
1545   pdgres[0]=211;pdgres[1]=2122;
1546   AliKFParticle *delta=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1547   Double_t probDe=TMath::Prob(delta->GetChi2(),delta->GetNDF());
1548   if(probDe>0.1) resNumber[2]=1;
1549
1550   return ;
1551
1552 }
1553 //-------------------------------------
1554 void AliAnalysisTaskSELambdac::IspKpiResonant(AliAODRecoDecayHF3Prong *d,Double_t field,Int_t *resNumber) const{
1555
1556   // apply PID using the resonant channels
1557   //if lambda* -> pk
1558   Double_t mass[2]={1.520,0.005};
1559   Int_t ipRes[2]={0,1};
1560   Int_t pdgres[2]={2212,321};
1561   AliKFParticle *lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1562   Double_t probLa=TMath::Prob(lambda1520->GetChi2(),lambda1520->GetNDF());
1563   if(probLa>0.1) resNumber[0]=1;
1564   //K* -> kpi
1565   mass[0]=0.8961;mass[1]=0.03;
1566   ipRes[0]=1;ipRes[1]=2;
1567   pdgres[1]=211;pdgres[0]=321;
1568   AliKFParticle *kstar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1569   Double_t probKa=TMath::Prob(kstar->GetChi2(),kstar->GetNDF());
1570   if(probKa>0.1) resNumber[1]=1;
1571
1572   //    Delta++
1573   mass[0]=1.232;mass[1]=0.15;
1574   ipRes[0]=0;ipRes[1]=2;
1575   pdgres[0]=2122;pdgres[1]=211;
1576   AliKFParticle *delta=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1577   Double_t probDe=TMath::Prob(delta->GetChi2(),delta->GetNDF());
1578   if(probDe>0.1) resNumber[2]=1;
1579
1580   return ;
1581
1582 }
1583 //---------------------------
1584 void AliAnalysisTaskSELambdac::FillMassHists(AliAODEvent *aod,AliAODRecoDecayHF3Prong *part,
1585                                                 TClonesArray *arrayMC, AliRDHFCutsLctopKpi *cuts,Int_t *nSelectedloose,Int_t *nSelectedtight)
1586 {
1587
1588   //if MC PID or no PID, unset PID
1589   if(!fRealPid) cuts->SetUsePID(kFALSE);
1590   Int_t selection=cuts->IsSelected(part,AliRDHFCuts::kCandidate,aod);
1591   if(selection>0){
1592     nSelectedloose[0]=nSelectedloose[0]+1;
1593     Int_t iPtBin = -1;
1594     Double_t ptCand = part->Pt();
1595     Int_t index=0;
1596
1597     for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
1598       if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
1599     }
1600
1601     Float_t pdgCode=-2;
1602     Float_t pdgCode1=-1;
1603     Float_t pdgCode2=-1;
1604     Int_t labDp=-1;
1605     Float_t deltaPx=0.;
1606     Float_t deltaPy=0.;
1607     Float_t deltaPz=0.;
1608     Float_t truePt=0.;
1609     Float_t xDecay=0.;
1610     Float_t yDecay=0.;
1611     Float_t zDecay=0.;
1612
1613     if(fReadMC){
1614       labDp = MatchToMCLambdac(part,arrayMC);
1615       if(labDp>0){
1616         AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
1617         AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
1618         AliAODMCParticle *dg1 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(1));
1619         deltaPx=partDp->Px()-part->Px();
1620         deltaPy=partDp->Py()-part->Py();
1621         deltaPz=partDp->Pz()-part->Pz();
1622         truePt=partDp->Pt();
1623         xDecay=dg0->Xv();
1624         yDecay=dg0->Yv();
1625         zDecay=dg0->Zv();
1626         pdgCode=TMath::Abs(partDp->GetPdgCode());
1627         pdgCode1=TMath::Abs(dg0->GetPdgCode());
1628         pdgCode2=TMath::Abs(dg1->GetPdgCode());
1629
1630       }else{
1631         pdgCode=-1;
1632       }
1633     }
1634
1635     Double_t invMasspKpi=-1.;
1636     Double_t invMasspiKp=-1.;
1637     Double_t invMasspKpiLpi=-1.;
1638     Double_t invMasspiKpLpi=-1.;
1639     Double_t invMasspKpiKp=-1.;
1640     Double_t invMasspiKpKp=-1.;
1641     Double_t invMasspiKpDk=-1.;
1642     Double_t invMasspKpiDk=-1.;
1643     Double_t invMasspiKp3Pr=-1.;
1644     Double_t invMasspKpi3Pr=-1.;
1645     Double_t field=aod->GetMagneticField();
1646     //apply MC PID
1647     if(fReadMC && fMCPid){
1648       if(IspKpiMC(part,arrayMC)) invMasspKpi=part->InvMassLcpKpi();
1649       if(IspiKpMC(part,arrayMC)) invMasspiKp=part->InvMassLcpiKp();
1650     }
1651
1652     // apply realistic PID
1653     if(fRealPid){
1654       if(selection==1 || selection==3) invMasspKpi=part->InvMassLcpKpi();
1655       if(selection>=2) invMasspiKp=part->InvMassLcpiKp();
1656     }
1657     //apply PID using resonances
1658     if (fResPid && fRealPid) {
1659       Int_t ispKpi[3]={0,0,0};
1660       IspKpiResonant(part,field,ispKpi);
1661       Int_t ispiKp[3]={0,0,0};
1662       IspiKpResonant(part,field,ispiKp);
1663       if (selection==3 || selection==1) {
1664        if(ispKpi[0]==1){
1665         invMasspKpiLpi=part->InvMassLcpKpi();
1666        }
1667        if(ispKpi[1]==1){
1668         invMasspKpiKp=part->InvMassLcpKpi();
1669        }
1670        if(ispKpi[2]==1){
1671         invMasspKpiDk=part->InvMassLcpKpi();
1672        }
1673        if(ispKpi[2]==0 && ispKpi[1]==0 && ispKpi[0]==0){
1674         invMasspKpi3Pr=part->InvMassLcpKpi();
1675        }
1676       }
1677       if(selection>=2) {
1678        if(ispiKp[0]==1){
1679         invMasspiKpLpi=part->InvMassLcpiKp();
1680        }
1681        if(ispiKp[1]==1){
1682         invMasspiKpKp=part->InvMassLcpiKp();
1683        }
1684        if(ispiKp[2]==1){
1685         invMasspiKpDk=part->InvMassLcpiKp();
1686        }
1687        if(ispiKp[2]==0 && ispiKp[1]==0 && ispiKp[0]==0){
1688         invMasspiKp3Pr=part->InvMassLcpiKp();
1689        }
1690       }
1691
1692      }
1693
1694
1695     if(invMasspiKp<0. && invMasspKpi<0.) return;
1696
1697     Int_t passTightCuts=0;
1698     if(fAnalysis) {
1699      passTightCuts=fRDCutsAnalysis->IsSelected(part,AliRDHFCuts::kCandidate,aod);
1700      if(fUseKF){
1701       Int_t pdgs[3]={0,0,0};
1702       if(invMasspKpi>0.){
1703        pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
1704        if(!VertexingKF(part,pdgs,field)) {
1705         invMasspKpi=-1.;
1706         invMasspKpi3Pr=-1.;
1707         invMasspKpiDk=-1.;
1708         invMasspKpiLpi=-1.;
1709         invMasspKpiKp=-1.;
1710        }
1711       }
1712       if(invMasspiKp>0.){
1713         pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
1714         if(!VertexingKF(part,pdgs,field)) {
1715           invMasspiKp=-1.;
1716           invMasspiKp3Pr=-1.;
1717           invMasspiKpDk=-1.;
1718           invMasspiKpLpi=-1.;
1719           invMasspiKpKp=-1.;
1720        }
1721      }
1722    }
1723     if(passTightCuts>0) nSelectedtight[0]=nSelectedtight[0]+1;
1724   }
1725
1726     Float_t tmp[24];
1727     if (fFillNtuple) {
1728       tmp[0]=pdgCode;
1729       tmp[1]=deltaPx;
1730       tmp[2]=deltaPy;
1731       tmp[3]=deltaPz;
1732       tmp[4]=truePt;
1733       tmp[5]=xDecay;
1734       tmp[6]=yDecay;
1735       tmp[7]=zDecay;
1736       if(pdgCode1==2212) {tmp[8]=part->PtProng(0);}else{tmp[8]=0.;}
1737       if(pdgCode1==211) {tmp[10]=part->PtProng(0);}else{tmp[10]=0.;}
1738       tmp[9]=part->PtProng(1);
1739       if(pdgCode2==211) {tmp[10]=part->PtProng(2);}else{tmp[10]=0.;}
1740       tmp[11]=part->Pt();
1741       tmp[12]=part->CosPointingAngle();
1742       tmp[13]=part->DecayLength();
1743       tmp[14]=part->Xv();
1744       tmp[15]=part->Yv();
1745       tmp[16]=part->Zv();
1746       if(invMasspiKp>0.) tmp[17]=invMasspiKp;
1747       if(invMasspKpi>0.) tmp[17]=invMasspKpi;
1748       tmp[18]=part->GetSigmaVert();
1749       tmp[19]=part->Getd0Prong(0);
1750       tmp[20]=part->Getd0Prong(1);
1751       tmp[21]=part->Getd0Prong(2);
1752       tmp[22]=part->GetDCA();
1753       tmp[23]=part->Prodd0d0();
1754       fNtupleLambdac->Fill(tmp);
1755       PostData(7,fNtupleLambdac);
1756     }
1757
1758     if(part->Pt()>3.&& part->Pt()<=6.){
1759       if(invMasspiKp>0. && invMasspKpi>0.){
1760         if(invMasspiKp>0.) fhMassPtGreater3->Fill(invMasspiKp,0.5);
1761         if(invMasspKpi>0.) fhMassPtGreater3->Fill(invMasspKpi,0.5);
1762       }else{
1763         if(invMasspiKp>0.) fhMassPtGreater3->Fill(invMasspiKp);
1764         if(invMasspKpi>0.) fhMassPtGreater3->Fill(invMasspKpi);
1765       }
1766       if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1767         if(invMasspiKpLpi>0.) fhMassPtGreater3Lpi->Fill(invMasspiKpLpi,0.5);
1768         if(invMasspKpiLpi>0.) fhMassPtGreater3Lpi->Fill(invMasspKpiLpi,0.5);
1769       }else{
1770         if(invMasspiKpLpi>0.) fhMassPtGreater3Lpi->Fill(invMasspiKpLpi);
1771         if(invMasspKpiLpi>0.) fhMassPtGreater3Lpi->Fill(invMasspKpiLpi);
1772       }
1773       if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1774         if(invMasspiKpKp>0.) fhMassPtGreater3Kp->Fill(invMasspiKpKp,0.5);
1775         if(invMasspKpiKp>0.) fhMassPtGreater3Kp->Fill(invMasspKpiKp,0.5);
1776       }else{
1777         if(invMasspiKpKp>0.) fhMassPtGreater3Kp->Fill(invMasspiKpKp);
1778         if(invMasspKpiKp>0.) fhMassPtGreater3Kp->Fill(invMasspKpiKp);
1779       }
1780       if(invMasspiKpDk>0. && invMasspKpiDk>0.){
1781        if(invMasspiKpDk>0.) fhMassPtGreater3Dk->Fill(invMasspiKpDk,0.5);
1782        if(invMasspKpiDk>0.) fhMassPtGreater3Dk->Fill(invMasspKpiDk,0.5);
1783       }else{
1784        if(invMasspiKpDk>0.) fhMassPtGreater3Dk->Fill(invMasspiKpDk);
1785        if(invMasspKpiDk>0.) fhMassPtGreater3Dk->Fill(invMasspKpiDk);
1786       }
1787       if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0.){
1788         if(invMasspiKp3Pr>0.) fhMassPtGreater33Pr->Fill(invMasspiKp3Pr,0.5);
1789         if(invMasspKpi3Pr>0.) fhMassPtGreater33Pr->Fill(invMasspKpi3Pr,0.5);
1790       }else{
1791         if(invMasspiKp3Pr>0.) fhMassPtGreater33Pr->Fill(invMasspiKp3Pr);
1792         if(invMasspKpi3Pr>0.) fhMassPtGreater33Pr->Fill(invMasspKpi3Pr);
1793       }
1794
1795       if(passTightCuts>0){
1796         if(invMasspiKp>0. && invMasspKpi>0.){
1797           if(invMasspiKp>0.) fhMassPtGreater3TC->Fill(invMasspiKp,0.5);
1798           if(invMasspKpi>0.) fhMassPtGreater3TC->Fill(invMasspKpi,0.5);
1799         }else{
1800           if(invMasspiKp>0.) fhMassPtGreater3TC->Fill(invMasspiKp);
1801           if(invMasspKpi>0.) fhMassPtGreater3TC->Fill(invMasspKpi);
1802         }
1803         if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1804           if(invMasspiKpLpi>0.) fhMassPtGreater3LpiTC->Fill(invMasspiKpLpi,0.5);
1805           if(invMasspKpiLpi>0.) fhMassPtGreater3LpiTC->Fill(invMasspKpiLpi,0.5);
1806         }else{
1807           if(invMasspiKpLpi>0.) fhMassPtGreater3LpiTC->Fill(invMasspiKpLpi);
1808           if(invMasspKpiLpi>0.) fhMassPtGreater3LpiTC->Fill(invMasspKpiLpi);
1809         }
1810         if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1811           if(invMasspiKpKp>0.) fhMassPtGreater3KpTC->Fill(invMasspiKpKp,0.5);
1812           if(invMasspKpiKp>0.) fhMassPtGreater3KpTC->Fill(invMasspKpiKp,0.5);
1813         }else{
1814           if(invMasspiKpKp>0.) fhMassPtGreater3KpTC->Fill(invMasspiKpKp);
1815           if(invMasspKpiKp>0.) fhMassPtGreater3KpTC->Fill(invMasspKpiKp);
1816         }
1817
1818        if(invMasspiKpDk>0. && invMasspKpiDk>0.){
1819         if(invMasspiKpDk>0.) fhMassPtGreater3DkTC->Fill(invMasspiKpDk,0.5);
1820         if(invMasspKpiDk>0.) fhMassPtGreater3DkTC->Fill(invMasspKpiDk,0.5);
1821        }else{
1822         if(invMasspiKpDk>0.) fhMassPtGreater3DkTC->Fill(invMasspiKpDk);
1823         if(invMasspKpiDk>0.) fhMassPtGreater3DkTC->Fill(invMasspKpiDk);
1824        }
1825       if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0.){
1826        if(invMasspiKp3Pr>0.) fhMassPtGreater33PrTC->Fill(invMasspiKp3Pr,0.5);
1827        if(invMasspKpi3Pr>0.) fhMassPtGreater33PrTC->Fill(invMasspKpi3Pr,0.5);
1828       }else{
1829        if(invMasspiKp3Pr>0.) fhMassPtGreater33PrTC->Fill(invMasspiKp3Pr);
1830        if(invMasspKpi3Pr>0.) fhMassPtGreater33PrTC->Fill(invMasspKpi3Pr);
1831       }
1832
1833      }
1834     }
1835     if(part->Pt()>2.&& part->Pt()<=6.){
1836       if(invMasspiKp>0. && invMasspKpi>0.){
1837         if(invMasspiKp>0.) fhMassPtGreater2->Fill(invMasspiKp,0.5);
1838         if(invMasspKpi>0.) fhMassPtGreater2->Fill(invMasspKpi,0.5);
1839       }else{
1840         if(invMasspiKp>0.) fhMassPtGreater2->Fill(invMasspiKp);
1841         if(invMasspKpi>0.) fhMassPtGreater2->Fill(invMasspKpi);
1842       }
1843       if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1844         if(invMasspiKpLpi>0.) fhMassPtGreater2Lpi->Fill(invMasspiKpLpi,0.5);
1845         if(invMasspKpiLpi>0.) fhMassPtGreater2Lpi->Fill(invMasspKpiLpi,0.5);
1846       }else{
1847         if(invMasspiKpLpi>0.) fhMassPtGreater2Lpi->Fill(invMasspiKpLpi);
1848         if(invMasspKpiLpi>0.) fhMassPtGreater2Lpi->Fill(invMasspKpiLpi);
1849       }
1850       if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1851         if(invMasspiKpKp>0.) fhMassPtGreater2Kp->Fill(invMasspiKpKp,0.5);
1852         if(invMasspKpiKp>0.) fhMassPtGreater2Kp->Fill(invMasspKpiKp,0.5);
1853       }else{
1854         if(invMasspiKpKp>0.) fhMassPtGreater2Kp->Fill(invMasspiKpKp);
1855         if(invMasspKpiKp>0.) fhMassPtGreater2Kp->Fill(invMasspKpiKp);
1856       }
1857       if(invMasspiKpDk>0. && invMasspKpiDk>0.){
1858        if(invMasspiKpDk>0.) fhMassPtGreater2Dk->Fill(invMasspiKpDk,0.5);
1859        if(invMasspKpiDk>0.) fhMassPtGreater2Dk->Fill(invMasspKpiDk,0.5);
1860       }else{
1861        if(invMasspiKpDk>0.) fhMassPtGreater2Dk->Fill(invMasspiKpDk);
1862        if(invMasspKpiDk>0.) fhMassPtGreater2Dk->Fill(invMasspKpiDk);
1863       }
1864      if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0.){
1865       if(invMasspiKp3Pr>0.) fhMassPtGreater23Pr->Fill(invMasspiKp3Pr,0.5);
1866       if(invMasspKpi3Pr>0.) fhMassPtGreater23Pr->Fill(invMasspKpi3Pr,0.5);
1867      }else{
1868       if(invMasspiKp3Pr>0.) fhMassPtGreater23Pr->Fill(invMasspiKp3Pr);
1869       if(invMasspKpi3Pr>0.) fhMassPtGreater23Pr->Fill(invMasspKpi3Pr);
1870     }
1871
1872       if(passTightCuts>0){
1873         if(invMasspiKp>0. && invMasspKpi>0.){
1874           if(invMasspiKp>0.) fhMassPtGreater2TC->Fill(invMasspiKp,0.5);
1875           if(invMasspKpi>0.) fhMassPtGreater2TC->Fill(invMasspKpi,0.5);
1876         }else{
1877           if(invMasspiKp>0.) fhMassPtGreater2TC->Fill(invMasspiKp);
1878           if(invMasspKpi>0.) fhMassPtGreater2TC->Fill(invMasspKpi);
1879         }
1880         if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1881           if(invMasspiKpLpi>0.) fhMassPtGreater2LpiTC->Fill(invMasspiKpLpi,0.5);
1882           if(invMasspKpiLpi>0.) fhMassPtGreater2LpiTC->Fill(invMasspKpiLpi,0.5);
1883         }else{
1884           if(invMasspiKpLpi>0.) fhMassPtGreater2LpiTC->Fill(invMasspiKpLpi);
1885           if(invMasspKpiLpi>0.) fhMassPtGreater2LpiTC->Fill(invMasspKpiLpi);
1886         }
1887         if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1888           if(invMasspiKpKp>0.) fhMassPtGreater2KpTC->Fill(invMasspiKpKp,0.5);
1889           if(invMasspKpiKp>0.) fhMassPtGreater2KpTC->Fill(invMasspKpiKp,0.5);
1890         }else{
1891           if(invMasspiKpKp>0.) fhMassPtGreater2KpTC->Fill(invMasspiKpKp);
1892           if(invMasspKpiKp>0.) fhMassPtGreater2KpTC->Fill(invMasspKpiKp);
1893         }
1894        if(invMasspiKpDk>0. && invMasspKpiDk>0.){
1895         if(invMasspiKpDk>0.) fhMassPtGreater2DkTC->Fill(invMasspiKpDk,0.5);
1896         if(invMasspKpiDk>0.) fhMassPtGreater2DkTC->Fill(invMasspKpiDk,0.5);
1897        }else{
1898         if(invMasspiKpDk>0.) fhMassPtGreater2DkTC->Fill(invMasspiKpDk);
1899         if(invMasspKpiDk>0.) fhMassPtGreater2DkTC->Fill(invMasspKpiDk);
1900        }
1901        if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0.){
1902         if(invMasspiKp3Pr>0.) fhMassPtGreater23PrTC->Fill(invMasspiKp3Pr,0.5);
1903         if(invMasspKpi3Pr>0.) fhMassPtGreater23PrTC->Fill(invMasspKpi3Pr,0.5);
1904       }else{
1905        if(invMasspiKp3Pr>0.) fhMassPtGreater23PrTC->Fill(invMasspiKp3Pr);
1906        if(invMasspKpi3Pr>0.) fhMassPtGreater23PrTC->Fill(invMasspKpi3Pr);
1907       }
1908
1909     }
1910    }
1911
1912     if(iPtBin>=0){
1913       index=GetHistoIndex(iPtBin);
1914       if(invMasspiKp>0. && invMasspKpi>0.){
1915         if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp,0.5);
1916         if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi,0.5);
1917       }else{
1918         if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
1919         if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
1920       }
1921       if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1922         if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi,0.5);
1923         if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi,0.5);
1924       }else{
1925         if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi);
1926         if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi);
1927       }
1928       if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1929         if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp,0.5);
1930         if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp,0.5);
1931       }else{
1932         if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp);
1933         if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp);
1934       }
1935       if(invMasspiKpDk>0. && invMasspKpiDk>0.){
1936        if(invMasspiKpDk>0.) fMassHistDk[index]->Fill(invMasspiKpDk,0.5);
1937        if(invMasspKpiDk>0.) fMassHistDk[index]->Fill(invMasspKpiDk,0.5);
1938       }else{
1939        if(invMasspiKpDk>0.) fMassHistDk[index]->Fill(invMasspiKpDk);
1940        if(invMasspKpiDk>0.) fMassHistDk[index]->Fill(invMasspKpiDk);
1941       }
1942       if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0.){
1943        if(invMasspiKp3Pr>0.) fMassHist3Pr[index]->Fill(invMasspiKp3Pr,0.5);
1944        if(invMasspKpi3Pr>0.) fMassHist3Pr[index]->Fill(invMasspKpi3Pr,0.5);
1945       }else{
1946        if(invMasspiKp3Pr>0.) fMassHist3Pr[index]->Fill(invMasspiKp3Pr);
1947        if(invMasspKpi3Pr>0.) fMassHist3Pr[index]->Fill(invMasspKpi3Pr);
1948      }
1949
1950       if(passTightCuts>0){
1951         if(invMasspiKp>0. && invMasspKpi>0. && passTightCuts==3){
1952           if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp,0.5);
1953           if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi,0.5);
1954         }else{
1955           if(invMasspiKp>0. && passTightCuts==2) fMassHistTC[index]->Fill(invMasspiKp);
1956           if(invMasspKpi>0. && passTightCuts==1) fMassHistTC[index]->Fill(invMasspKpi);
1957         }
1958         if(invMasspiKpLpi>0. && invMasspKpiLpi>0. && passTightCuts==3){
1959           if(invMasspiKpLpi>0.) fMassHistLpiTC[index]->Fill(invMasspiKpLpi,0.5);
1960           if(invMasspKpiLpi>0.) fMassHistLpiTC[index]->Fill(invMasspKpiLpi,0.5);
1961         }else{
1962           if(invMasspiKpLpi>0. && passTightCuts==2) fMassHistLpiTC[index]->Fill(invMasspiKpLpi);
1963           if(invMasspKpiLpi>0.&& passTightCuts==1) fMassHistLpiTC[index]->Fill(invMasspKpiLpi);
1964         }
1965         if(invMasspiKpKp>0. && invMasspKpiKp>0. && passTightCuts==3){
1966           if(invMasspiKpKp>0.) fMassHistKpTC[index]->Fill(invMasspiKpKp,0.5);
1967           if(invMasspKpiKp>0.) fMassHistKpTC[index]->Fill(invMasspKpiKp,0.5);
1968         }else{
1969           if(invMasspiKpKp>0. && passTightCuts==2) fMassHistKpTC[index]->Fill(invMasspiKpKp);
1970           if(invMasspKpiKp>0.&& passTightCuts==1) fMassHistKpTC[index]->Fill(invMasspKpiKp);
1971         }
1972        if(invMasspiKpDk>0. && invMasspKpiDk>0. && passTightCuts==3){
1973         if(invMasspiKpDk>0.) fMassHistDkTC[index]->Fill(invMasspiKpDk,0.5);
1974         if(invMasspKpiDk>0.) fMassHistDkTC[index]->Fill(invMasspKpiDk,0.5);
1975        }else{
1976         if(invMasspiKpDk>0. && passTightCuts==2) fMassHistDkTC[index]->Fill(invMasspiKpDk);
1977         if(invMasspKpiDk>0.&& passTightCuts==1) fMassHistDkTC[index]->Fill(invMasspKpiDk);
1978        }
1979        if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0. && passTightCuts==3){
1980         if(invMasspiKp3Pr>0.) fMassHist3PrTC[index]->Fill(invMasspiKp3Pr,0.5);
1981         if(invMasspKpi3Pr>0.) fMassHist3PrTC[index]->Fill(invMasspKpi3Pr,0.5);
1982        }else{
1983         if(invMasspiKp3Pr>0. && passTightCuts==2) fMassHist3PrTC[index]->Fill(invMasspiKp3Pr);
1984         if(invMasspKpi3Pr>0.&& passTightCuts==1) fMassHist3PrTC[index]->Fill(invMasspKpi3Pr);
1985       }
1986
1987       }
1988
1989       if(fReadMC){
1990         if(labDp>0) {
1991           index=GetSignalHistoIndex(iPtBin);
1992           if(invMasspiKp>0. && invMasspKpi>0.){
1993             if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp,0.5);
1994             if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi,0.5);
1995           }else{
1996             if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
1997             if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
1998           }
1999           if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
2000             if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi,0.5);
2001             if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi,0.5);
2002           }else{
2003             if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi);
2004             if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi);
2005           }
2006           if(invMasspiKpKp>0. && invMasspKpiKp>0.){
2007             if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp,0.5);
2008             if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp,0.5);
2009           }else{
2010             if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp);
2011             if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp);
2012           }
2013           if(invMasspiKpDk>0. && invMasspKpiDk>0.){
2014             if(invMasspiKpDk>0.) fMassHistDk[index]->Fill(invMasspiKpDk,0.5);
2015             if(invMasspKpiDk>0.) fMassHistDk[index]->Fill(invMasspKpiDk,0.5);
2016           }else{
2017             if(invMasspiKpDk>0.) fMassHistDk[index]->Fill(invMasspiKpDk);
2018             if(invMasspKpiDk>0.) fMassHistDk[index]->Fill(invMasspKpiDk);
2019           }
2020           if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0.){
2021             if(invMasspiKp3Pr>0.) fMassHist3Pr[index]->Fill(invMasspiKp3Pr,0.5);
2022             if(invMasspKpi3Pr>0.) fMassHist3Pr[index]->Fill(invMasspKpi3Pr,0.5);
2023           }else{
2024             if(invMasspiKp3Pr>0.) fMassHist3Pr[index]->Fill(invMasspiKp3Pr);
2025             if(invMasspKpi3Pr>0.) fMassHistDk[index]->Fill(invMasspKpi3Pr);
2026           }
2027
2028           if(passTightCuts>0){
2029             if(invMasspiKp>0. && invMasspKpi>0. && passTightCuts==3){
2030               if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp,0.5);
2031               if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi,0.5);
2032             }else{
2033               if(invMasspiKp>0. && passTightCuts==2) fMassHistTC[index]->Fill(invMasspiKp);
2034               if(invMasspKpi>0.&& passTightCuts==1) fMassHistTC[index]->Fill(invMasspKpi);
2035             }
2036             if(invMasspiKpLpi>0. && invMasspKpiLpi>0. && passTightCuts==3){
2037               if(invMasspiKpLpi>0.) fMassHistLpiTC[index]->Fill(invMasspiKpLpi,0.5);
2038               if(invMasspKpiLpi>0.) fMassHistLpiTC[index]->Fill(invMasspKpiLpi,0.5);
2039             }else{
2040               if(invMasspiKpLpi>0. && passTightCuts==2) fMassHistLpiTC[index]->Fill(invMasspiKpLpi);
2041               if(invMasspKpiLpi>0.&& passTightCuts==1) fMassHistLpiTC[index]->Fill(invMasspKpiLpi);
2042             }
2043             if(invMasspiKpKp>0. && invMasspKpiKp>0. && passTightCuts==3){
2044               if(invMasspiKpKp>0.) fMassHistKpTC[index]->Fill(invMasspiKpKp,0.5);
2045               if(invMasspKpiKp>0.) fMassHistKpTC[index]->Fill(invMasspKpiKp,0.5);
2046             }else{
2047               if(invMasspiKpKp>0. && passTightCuts==2) fMassHistKpTC[index]->Fill(invMasspiKpKp);
2048               if(invMasspKpiKp>0.&& passTightCuts==1) fMassHistKpTC[index]->Fill(invMasspKpiKp);
2049             }
2050             if(invMasspiKpDk>0. && invMasspKpiDk>0. && passTightCuts==3){
2051               if(invMasspiKpDk>0.) fMassHistDkTC[index]->Fill(invMasspiKpDk,0.5);
2052               if(invMasspKpiDk>0.) fMassHistDkTC[index]->Fill(invMasspKpiDk,0.5);
2053             }else{
2054               if(invMasspiKpDk>0. && passTightCuts==2) fMassHistDkTC[index]->Fill(invMasspiKpDk);
2055               if(invMasspKpiDk>0.&& passTightCuts==1) fMassHistDkTC[index]->Fill(invMasspKpiDk);
2056             }
2057             if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0. && passTightCuts==3){
2058               if(invMasspiKp3Pr>0.) fMassHist3PrTC[index]->Fill(invMasspiKp3Pr,0.5);
2059               if(invMasspKpi3Pr>0.) fMassHist3PrTC[index]->Fill(invMasspKpi3Pr,0.5);
2060             }else{
2061               if(invMasspiKp3Pr>0. && passTightCuts==2) fMassHist3PrTC[index]->Fill(invMasspiKp3Pr);
2062               if(invMasspKpi3Pr>0.&& passTightCuts==1) fMassHist3PrTC[index]->Fill(invMasspKpi3Pr);
2063             }
2064           }
2065
2066         }else{
2067           index=GetBackgroundHistoIndex(iPtBin);
2068           if(invMasspiKp>0. && invMasspKpi>0.){
2069             fMassHist[index]->Fill(invMasspiKp,0.5);
2070             fMassHist[index]->Fill(invMasspKpi,0.5);
2071           }else{
2072             if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
2073             if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
2074           }
2075           if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
2076             if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi,0.5);
2077             if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi,0.5);
2078           }else{
2079             if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi);
2080             if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi);
2081           }
2082           if(invMasspiKpKp>0. && invMasspKpiKp>0.){
2083             if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp,0.5);
2084             if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp,0.5);
2085           }else{
2086             if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp);
2087             if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp);
2088           }
2089           if(invMasspiKpDk>0. && invMasspKpiDk>0.){
2090             if(invMasspiKpDk>0.) fMassHistDk[index]->Fill(invMasspiKpDk,0.5);
2091             if(invMasspKpiDk>0.) fMassHistDk[index]->Fill(invMasspKpiDk,0.5);
2092           }else{
2093             if(invMasspiKpDk>0.) fMassHistDk[index]->Fill(invMasspiKpDk);
2094             if(invMasspKpiDk>0.) fMassHistDk[index]->Fill(invMasspKpiDk);
2095           }
2096           if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0.){
2097             if(invMasspiKp3Pr>0.) fMassHist3Pr[index]->Fill(invMasspiKp3Pr,0.5);
2098             if(invMasspKpi3Pr>0.) fMassHist3Pr[index]->Fill(invMasspKpi3Pr,0.5);
2099           }else{
2100             if(invMasspiKp3Pr>0.) fMassHist3Pr[index]->Fill(invMasspiKp3Pr);
2101             if(invMasspKpi3Pr>0.) fMassHistDk[index]->Fill(invMasspKpi3Pr);
2102           }
2103           if(invMasspiKp>0. && invMasspKpi>0. && passTightCuts==3){
2104             fMassHistTC[index]->Fill(invMasspiKp,0.5);
2105             fMassHistTC[index]->Fill(invMasspKpi,0.5);
2106           }else{
2107             if(invMasspiKp>0. && passTightCuts==2) fMassHistTC[index]->Fill(invMasspiKp);
2108             if(invMasspKpi>0. && passTightCuts==1) fMassHistTC[index]->Fill(invMasspKpi);
2109           }
2110           if(invMasspiKpLpi>0. && invMasspKpiLpi>0. && passTightCuts==3){
2111             if(invMasspiKpLpi>0.) fMassHistLpiTC[index]->Fill(invMasspiKpLpi,0.5);
2112             if(invMasspKpiLpi>0.) fMassHistLpiTC[index]->Fill(invMasspKpiLpi,0.5);
2113           }else{
2114             if(invMasspiKpLpi>0. && passTightCuts==2) fMassHistLpiTC[index]->Fill(invMasspiKpLpi);
2115             if(invMasspKpiLpi>0.&& passTightCuts==1) fMassHistLpiTC[index]->Fill(invMasspKpiLpi);
2116           }
2117           if(invMasspiKpKp>0. && invMasspKpiKp>0. && passTightCuts==3){
2118             if(invMasspiKpKp>0.) fMassHistKpTC[index]->Fill(invMasspiKpKp,0.5);
2119             if(invMasspKpiKp>0.) fMassHistKpTC[index]->Fill(invMasspKpiKp,0.5);
2120           }else{
2121             if(invMasspiKpKp>0. && passTightCuts==2) fMassHistKpTC[index]->Fill(invMasspiKpKp);
2122             if(invMasspKpiKp>0.&& passTightCuts==1) fMassHistKpTC[index]->Fill(invMasspKpiKp);
2123           }
2124           if(invMasspiKpDk>0. && invMasspKpiDk>0. && passTightCuts==3){
2125             if(invMasspiKpDk>0.) fMassHistDkTC[index]->Fill(invMasspiKpDk,0.5);
2126             if(invMasspKpiDk>0.) fMassHistDkTC[index]->Fill(invMasspKpiDk,0.5);
2127           }else{
2128             if(invMasspiKpDk>0. && passTightCuts==2) fMassHistDkTC[index]->Fill(invMasspiKpDk);
2129             if(invMasspKpiDk>0.&& passTightCuts==1) fMassHistDkTC[index]->Fill(invMasspKpiDk);
2130            }
2131            if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0. && passTightCuts==3){
2132             if(invMasspiKp3Pr>0.) fMassHist3PrTC[index]->Fill(invMasspiKp3Pr,0.5);
2133             if(invMasspKpi3Pr>0.) fMassHist3PrTC[index]->Fill(invMasspKpi3Pr,0.5);
2134            }else{
2135             if(invMasspiKp3Pr>0. && passTightCuts==2) fMassHist3PrTC[index]->Fill(invMasspiKp3Pr);
2136             if(invMasspKpi3Pr>0.&& passTightCuts==1) fMassHist3PrTC[index]->Fill(invMasspKpi3Pr);
2137           }
2138         }
2139
2140       }
2141     }
2142   }
2143   return;
2144 }
2145 //-----------------------
2146 void AliAnalysisTaskSELambdac::FillVarHists(AliAODRecoDecayHF3Prong *part,
2147                                                TClonesArray *arrMC,
2148                                                AliRDHFCutsLctopKpi *cuts,
2149                                                /* TList *listout,*/
2150                                                AliAODEvent* aod) {
2151   //
2152   // function used in UserExec to fill variable histograms analysing MC
2153   //
2154
2155   TString fillthis="";
2156
2157   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
2158   Double_t invmasscut=0.05;
2159
2160   Int_t pdgDgLctopKpi[3]={2212,321,211};
2161   Int_t lab=-9999;
2162   if(fReadMC)
2163     lab=part->MatchToMC(4122,arrMC,3,pdgDgLctopKpi); //return MC particle label if the array corresponds to a Lc, -1 if not (cf. AliAODRecoDecay.cxx)
2164
2165   Int_t isSelectedPID=cuts->IsSelectedPID(part); // 0 rejected, 1 Lc -> p K- pi+ (K at center because different sign is there),
2166                                                  // 2 Lc -> pi+ K- p (K at center because different sign is there),
2167                                                  // 3 both (it should never happen...)
2168
2169   if (isSelectedPID==0)fNentries->Fill(7);
2170   if (isSelectedPID==1)fNentries->Fill(8);
2171   if (isSelectedPID==2)fNentries->Fill(9);
2172   if (isSelectedPID==3)fNentries->Fill(10);
2173
2174   AliAODTrack *prong0=(AliAODTrack*)part->GetDaughter(0);
2175   AliAODTrack *prong1=(AliAODTrack*)part->GetDaughter(1);
2176   AliAODTrack *prong2=(AliAODTrack*)part->GetDaughter(2);
2177   if (!prong0 || !prong1 || !prong2) {
2178     fNentries->Fill(5);
2179     return;
2180   }
2181
2182   Double_t minvLcpKpi = part->InvMassLcpKpi();
2183   Double_t minvLcpiKp = part->InvMassLcpiKp();
2184
2185   //check pdg of the prongs
2186   Int_t labprong[3]={-1,-1,-1};
2187   if(fReadMC){
2188     labprong[0]=prong0->GetLabel();
2189     labprong[1]=prong1->GetLabel();
2190     labprong[2]=prong2->GetLabel();
2191   }
2192
2193   AliAODMCParticle *mcprong=0;
2194   Int_t pdgProngMC[3]={-1,-1,-1};
2195   if(fReadMC) {
2196     for (Int_t iprong=0;iprong<3;iprong++){
2197       if(labprong[iprong]<0) continue;
2198       mcprong = (AliAODMCParticle*)arrMC->At(labprong[iprong]);
2199       pdgProngMC[iprong]=TMath::Abs(mcprong->GetPdgCode());
2200     }
2201   }
2202
2203   Int_t pdgProngPID[3]={-1,-1,-1};
2204   if(isSelectedPID>0){
2205     pdgProngPID[1]=321;
2206     if(isSelectedPID==1) {pdgProngPID[0]=2212;pdgProngPID[2]=211;}
2207     if(isSelectedPID==2) {pdgProngPID[0]=211;pdgProngPID[2]=2212;}
2208   }
2209
2210   //fill hRealTot ---> PID efficiency denominators
2211   cuts->SetUsePID(kFALSE); //PAOLA
2212   Int_t selectionCandidate=cuts->IsSelected(part,AliRDHFCuts::kCandidate,aod);//PAOLA
2213   if(fReadMC && selectionCandidate>0) { // 3prongs candidate x Lc (no PID)
2214     if(lab>0){ // Signal
2215       for (Int_t iprong=0; iprong<3; iprong++) {
2216         switch (pdgProngMC[iprong]) {
2217         case 2212:
2218           fillthis="hpRealTot";
2219           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2220           break;
2221         case 321:
2222           fillthis="hKRealTot";
2223           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2224           break;
2225         case 211:
2226           fillthis="hpiRealTot";
2227           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2228           break;
2229         default:
2230           break;
2231         }
2232       }
2233     } else { // bkg
2234       for (Int_t iprong=0; iprong<3; iprong++) {
2235         switch (pdgProngMC[iprong]) {
2236         case 2212:
2237           fillthis="hbpRealTot";
2238           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2239           break;
2240         case 321:
2241           fillthis="hbKRealTot";
2242           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2243           break;
2244         case 211:
2245           fillthis="hbpiRealTot";
2246           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2247           break;
2248         default:
2249           break;
2250         }
2251       }
2252     }
2253   }
2254
2255
2256   cuts->SetUsePID(kTRUE); //PAOLA
2257   Int_t selection=cuts->IsSelected(part,AliRDHFCuts::kCandidate,aod);
2258
2259   if ( (lab>0 && fReadMC) ||             // signal X MC
2260        (isSelectedPID>0 && !fReadMC) ) { // signal+bkg X real data
2261
2262     fillthis="hMass";
2263     if ( (fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 4122) ||
2264          (!fReadMC && (isSelectedPID>0 && part->Charge()>0)) ) { // Lc
2265       if ( (pdgProngPID[0]==2212) && (pdgProngPID[1]==321) && (pdgProngPID[2]==211) )
2266         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(minvLcpKpi); // signal+bkg X MC and real data
2267       else if ( (pdgProngPID[0]==211) && (pdgProngPID[1]==321) && (pdgProngPID[2]==2212) )
2268         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(minvLcpiKp); // signal+bkg X MC and real data
2269     }
2270     else if ( (fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == -4122) ||
2271               (!fReadMC && (isSelectedPID>0 && part->Charge()<0)) ) { // anti-Lc
2272       if ( (pdgProngPID[0]==2212) && (pdgProngPID[1]==321) && (pdgProngPID[2]==211) )
2273         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(minvLcpKpi); // signal+bkg X MC and real data
2274       else if ( (pdgProngPID[0]==211) && (pdgProngPID[1]==321) && (pdgProngPID[2]==2212) )
2275         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(minvLcpiKp); // signal+bkg X MC and real data
2276     }
2277
2278     if (selection>0) { // 3prongs candidate x Lc (yes PID)
2279
2280       Double_t ptmax=0.;
2281
2282       for (Int_t iprong=0; iprong<3; iprong++) {
2283         if (part->PtProng(iprong)>ptmax) ptmax=part->PtProng(iprong);
2284
2285         AliAODTrack *prong = (AliAODTrack*)part->GetDaughter(iprong);
2286         AliAODPid *pidObjtrk = (AliAODPid*)prong->GetDetPid();
2287         AliAODPidHF *pidObj = (AliAODPidHF*)cuts->GetPidHF();
2288         Bool_t hasTOF=pidObj->CheckStatus(prong,"TOF");
2289         Bool_t hasTPC=pidObj->CheckStatus(prong,"TPC");
2290         Double_t tofSignal=0.;
2291         Double_t dedxTPC=0.;
2292         Double_t momTOF=0.;
2293         Double_t momTPC=0.;
2294         if(hasTOF) {
2295           momTOF = prong->P();
2296           tofSignal=pidObjtrk->GetTOFsignal();
2297         }
2298         if(hasTPC) {
2299           momTPC = pidObjtrk->GetTPCmomentum();
2300           dedxTPC=pidObjtrk->GetTPCsignal();
2301         }
2302         switch (pdgProngPID[iprong]) {
2303         case 2212:
2304           fillthis="hpTOFSignal";
2305           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(tofSignal);
2306           fillthis="hpTPCSignal";
2307           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(dedxTPC);
2308           fillthis="hpptProng";
2309           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2310           fillthis="hpd0Prong";
2311           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
2312           fillthis="hpSignalVspTPC";
2313           ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,dedxTPC);
2314           fillthis="hpSignalVspTOF";
2315           ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,tofSignal);
2316           AliPID::EParticleType typep;
2317           typep=AliPID::EParticleType(4);
2318           if(hasTPC) {
2319             //Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typep);
2320             Double_t nsigmap = fPIDResponse->NumberOfSigmasTPC(prong,typep);
2321             fillthis="hpSigmaVspTPC";
2322             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,nsigmap);
2323           }
2324           if(hasTOF){
2325             //Double_t nsigma=fUtilPid->NumberOfSigmasTOF(prong,typep);
2326             Double_t nsigma=fPIDResponse->NumberOfSigmasTOF(prong,typep);
2327             fillthis="hpSigmaVspTOF";
2328             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,nsigma);
2329           }
2330
2331           if (fReadMC) { // fill hpIDTot ---> PID contamination denominator
2332                          //      hIDGood, hnoID ---> PID numerators
2333             fillthis="hpIDTot";
2334             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2335             if(pdgProngMC[iprong]==2212) { // id protons
2336               fillthis="hpIDGood";
2337               ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2338             }
2339             else { // misidentified electrons, muons, pions and kaons
2340               fillthis="hnopIDp";
2341               ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2342             }
2343           }
2344           break;
2345
2346         case 321:
2347           fillthis="hKTOFSignal";
2348           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(tofSignal);
2349           fillthis="hKTPCSignal";
2350           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(dedxTPC);
2351           fillthis="hKptProng";
2352           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2353           fillthis="hKd0Prong";
2354           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
2355           fillthis="hKSignalVspTPC";
2356           ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,dedxTPC);
2357           fillthis="hKSignalVspTOF";
2358           ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,tofSignal);
2359           AliPID::EParticleType typek;
2360           typek=AliPID::EParticleType(3);
2361           if(hasTPC) {
2362             //Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typek);
2363             Double_t nsigmap = fPIDResponse->NumberOfSigmasTPC(prong,typek);
2364             fillthis="hKSigmaVspTPC";
2365             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,nsigmap);
2366           }
2367           if(hasTOF){
2368             //Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typek);
2369             Double_t nsigma=fPIDResponse->NumberOfSigmasTOF(prong,typek); // protection against 'old' AODs...
2370             fillthis="hKSigmaVspTOF";
2371             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,nsigma);
2372           }
2373
2374           if (fReadMC) { // fill hKIDTot ---> PID contamination denominator
2375                          //      hIDGood, hnoID ---> PID numerators
2376             fillthis="hKIDTot";
2377             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2378             if(pdgProngMC[iprong]==321) { // id kaons
2379               fillthis="hKIDGood";
2380               ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2381             }
2382             else { // misidentified electrons, muons, pions and protons
2383               fillthis="hnokIDk";
2384               ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2385             }
2386           }
2387           break;
2388
2389         case 211:
2390           fillthis="hpiTOFSignal";
2391           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(tofSignal);
2392           fillthis="hpiTPCSignal";
2393           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(dedxTPC);
2394           fillthis="hpiptProng";
2395           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2396           fillthis="hpid0Prong";
2397           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
2398           fillthis="hpiSignalVspTPC";
2399           ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,dedxTPC);
2400           fillthis="hpiSignalVspTOF";
2401           ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,tofSignal);
2402           AliPID::EParticleType typepi;
2403           typepi=AliPID::EParticleType(2);
2404           if(hasTPC) {
2405             //Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typepi);
2406             Double_t nsigmap = fPIDResponse->NumberOfSigmasTPC(prong,typepi);
2407             fillthis="hpiSigmaVspTPC";
2408             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,nsigmap);
2409           }
2410           if(hasTOF){
2411             //Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typepi);
2412             Double_t nsigma=fPIDResponse->NumberOfSigmasTOF(prong,typepi); // protection against 'old' AODs...
2413             fillthis="hpiSigmaVspTOF";
2414             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,nsigma);
2415           }
2416
2417           if (fReadMC) { // fill hpiIDTot ---> PID contamination denominator
2418                          //      hIDGood, hnoID ---> PID numerators
2419             fillthis="hpiIDTot";
2420             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2421             if(pdgProngMC[iprong]==211) { // id pions
2422               fillthis="hpiIDGood";
2423               ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2424             }
2425             else { // misidentified electrons, muons, kaons and protons
2426               fillthis="hnopiIDpi";
2427               ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2428             }
2429           }
2430           break;
2431
2432         default:
2433           break;
2434         }
2435
2436       } //end loop on prongs
2437
2438       //now histograms where we don't need to check identity
2439       fillthis = "hDist12toPrim";
2440       ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->GetDist12toPrim());
2441       ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->GetDist23toPrim());
2442       fillthis = "hSigmaVert";
2443       ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->GetSigmaVert());
2444       fillthis = "hDCAs";
2445       Double_t dcas[3];
2446       part->GetDCAs(dcas);
2447       for (Int_t idca=0;idca<3;idca++) ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(dcas[idca]);
2448       fillthis = "hCosPointingAngle";
2449       ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->CosPointingAngle());
2450       fillthis = "hDecayLength";
2451       ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->DecayLength());
2452       Double_t sum2=part->Getd0Prong(0)*part->Getd0Prong(0)+
2453         part->Getd0Prong(1)*part->Getd0Prong(1)+
2454         part->Getd0Prong(2)*part->Getd0Prong(2);
2455       fillthis = "hSum2";
2456       ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(sum2);
2457       fillthis = "hptmax";
2458       ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(ptmax);
2459       fillthis="hLcpt";
2460       ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->Pt());
2461
2462     } // end if (selection)
2463
2464   } else if( lab<=0 && fReadMC) { // bkg x MC
2465
2466     fillthis="hbMass";
2467
2468     if (part->Charge()>0) {      //Lc background
2469       if ( (pdgProngPID[0]==2212) && (pdgProngPID[1]==321) && (pdgProngPID[2]==211) )
2470         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(minvLcpKpi); // bkg X MC
2471       else if ( (pdgProngPID[0]==211) && (pdgProngPID[1]==321) && (pdgProngPID[2]==2212) )
2472         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(minvLcpiKp); // bkg X MC
2473     }
2474     else if (part->Charge()<0){ //anti-Lc background
2475       if ( (pdgProngPID[0]==2212) && (pdgProngPID[1]==321) && (pdgProngPID[2]==211) )
2476         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(minvLcpKpi); // bkg X MC
2477       else if ( (pdgProngPID[0]==211) && (pdgProngPID[1]==321) && (pdgProngPID[2]==2212) )
2478         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(minvLcpiKp); // bkg X MC
2479     }
2480
2481
2482     //apply cut on invariant mass on the pair
2483     if (TMath::Abs(minvLcpKpi-mPDG)<invmasscut || TMath::Abs(minvLcpiKp-mPDG)<invmasscut) {
2484
2485       if (selection>0) { // 3prongs candidate x Lc (yes PID)
2486
2487         Double_t ptmax=0.;
2488         for (Int_t iprong=0; iprong<3; iprong++) {
2489           if(part->PtProng(iprong)>ptmax)ptmax=part->PtProng(iprong);
2490
2491           AliAODTrack *prong = (AliAODTrack*)part->GetDaughter(iprong);
2492           AliAODPid *pidObjtrk = (AliAODPid*)prong->GetDetPid();
2493           AliAODPidHF *pidObj = (AliAODPidHF*)cuts->GetPidHF();
2494           Bool_t hasTOF=pidObj->CheckStatus(prong,"TOF");
2495           Bool_t hasTPC=pidObj->CheckStatus(prong,"TPC");
2496           Double_t tofSignal=0.;
2497           Double_t dedxTPC=0.;
2498           Double_t momTOF=0.;
2499           Double_t momTPC=0.;
2500           if(hasTOF) {
2501             momTOF = prong->P();
2502             tofSignal=pidObjtrk->GetTOFsignal();
2503           }
2504           if(hasTPC) {
2505             momTPC = pidObjtrk->GetTPCmomentum();
2506             dedxTPC=pidObjtrk->GetTPCsignal();
2507           }
2508
2509           switch (pdgProngPID[iprong]) {
2510           case 2212:
2511             fillthis="hbpTOFSignal";
2512             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(tofSignal);
2513             fillthis="hbpTPCSignal";
2514             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(dedxTPC);
2515             fillthis="hbpptProng";
2516             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2517             fillthis="hbpd0Prong";
2518             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
2519             fillthis="hbpSignalVspTPC";
2520             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,dedxTPC);
2521             fillthis="hbpSignalVspTOF";
2522             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,tofSignal);
2523             AliPID::EParticleType typep;
2524             typep=AliPID::EParticleType(4);
2525             if(hasTPC) {
2526               //Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typep);
2527               Double_t nsigmap = fPIDResponse->NumberOfSigmasTPC(prong,typep);
2528               fillthis="hbpSigmaVspTPC";
2529               ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,nsigmap);
2530             }
2531             if(hasTOF){
2532               //Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typep);
2533               Double_t nsigma=fPIDResponse->NumberOfSigmasTOF(prong,typep); // protection against 'old' AODs...
2534               fillthis="hbpSigmaVspTOF";
2535               ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,nsigma);
2536             }
2537             if(fReadMC){ // fill hbpIDTot ---> PID contamination denominator
2538                          //      hbIDGood, hbnoID ---> PID numerators
2539               fillthis="hbpIDTot";
2540               ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2541               if(pdgProngMC[iprong]==2212) { // id protons
2542                 fillthis="hbpIDGood";
2543                 ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2544               }
2545               else { // misidentified electrons, muons, pions and kaons
2546                 fillthis="hbnopIDp";
2547                 ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2548               }
2549             }
2550             break;
2551
2552           case 321:
2553             fillthis="hbKTOFSignal";
2554             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(tofSignal);
2555             fillthis="hbKTPCSignal";
2556             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(dedxTPC);
2557             fillthis="hbKptProng";
2558             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2559             fillthis="hbKd0Prong";
2560             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
2561             fillthis="hbKSignalVspTPC";
2562             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,dedxTPC);
2563             fillthis="hbKSignalVspTOF";
2564             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,tofSignal);
2565             AliPID::EParticleType typek;
2566             typek=AliPID::EParticleType(3);
2567             if(hasTPC) {
2568               //Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typek);
2569               Double_t nsigmap = fPIDResponse->NumberOfSigmasTPC(prong,typek);
2570               fillthis="hbKSigmaVspTPC";
2571               ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,nsigmap);
2572             }
2573             if(hasTOF){
2574               //Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typek);
2575               Double_t nsigma=fPIDResponse->NumberOfSigmasTOF(prong,typek); // protection against 'old' AODs...
2576               fillthis="hbKSigmaVspTOF";
2577               ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,nsigma);
2578             }
2579             if (fReadMC) { // fill hbKIDTot ---> PID contamination denominator
2580                            //      hIDGood, hnoID ---> PID numerators
2581               fillthis="hbKIDTot";
2582               ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2583               if(pdgProngMC[iprong]==321) { // id kaons
2584                 fillthis="hbKIDGood";
2585                 ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2586               }
2587               else { // misidentified electrons, muons, pions and protons
2588                 fillthis="hbnokIDk";
2589                 ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2590               }
2591             }
2592             break;
2593
2594           case 211:
2595             fillthis="hbpiTOFSignal";
2596             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(tofSignal);
2597             fillthis="hbpiTPCSignal";
2598             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(dedxTPC);
2599             fillthis="hbpiptProng";
2600             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2601             fillthis="hbpid0Prong";
2602             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
2603             fillthis="hbpiSignalVspTPC";
2604             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,dedxTPC);
2605             fillthis="hbpiSignalVspTOF";
2606             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,tofSignal);
2607             AliPID::EParticleType typepi;
2608             typepi=AliPID::EParticleType(2);
2609             if(hasTPC) {
2610               //Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typepi);
2611               Double_t nsigmap = fPIDResponse->NumberOfSigmasTPC(prong,typepi);
2612               fillthis="hbpiSigmaVspTPC";
2613               ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,nsigmap);
2614             }
2615             if(hasTOF){
2616               //Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typepi);
2617               Double_t nsigma=fPIDResponse->NumberOfSigmasTOF(prong,typepi); // protection against 'old' AODs...
2618               fillthis="hbpiSigmaVspTOF";
2619               ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,nsigma);
2620             }
2621             if (fReadMC) { // fill hbpiIDTot ---> PID contamination denominator
2622                            //      hIDGood, hnonIDTot ---> PID numerators
2623               fillthis="hbpiIDTot";
2624               ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2625               if(pdgProngMC[iprong]==211) { // id pions
2626                 fillthis="hbpiIDGood";
2627                 ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2628               }
2629               else { // misidentified electrons, muons, kaons and protons
2630                 fillthis="hbnopiIDpi";
2631                 ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2632               }
2633             }
2634             break;
2635
2636           default:
2637             break;
2638           }
2639
2640         } //end loop on prongs
2641
2642         //now histograms where we don't need to check identity
2643         fillthis="hbLcpt";
2644         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->Pt());
2645
2646         fillthis = "hbDist12toPrim";
2647         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->GetDist12toPrim());
2648         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->GetDist23toPrim());
2649         fillthis = "hbSigmaVert";
2650         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->GetSigmaVert());
2651         fillthis = "hbDCAs";
2652         Double_t dcas[3];
2653         part->GetDCAs(dcas);
2654         for (Int_t idca=0;idca<3;idca++)
2655           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(dcas[idca]);
2656         fillthis = "hbCosPointingAngle";
2657         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->CosPointingAngle());
2658         fillthis = "hbDecayLength";
2659         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->DecayLength());
2660         Double_t sum2=part->Getd0Prong(0)*part->Getd0Prong(0)+
2661           part->Getd0Prong(1)*part->Getd0Prong(1)+
2662           part->Getd0Prong(2)*part->Getd0Prong(2);
2663         fillthis = "hbSum2";
2664         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(sum2);
2665         fillthis = "hbptmax";
2666         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(ptmax);
2667
2668       } // end if (selection)
2669
2670     } // end mass cut
2671
2672   } // end background case
2673
2674   return;
2675
2676 }
2677
2678 //-------------------------
2679 void AliAnalysisTaskSELambdac::FillAPrioriConcentrations(AliAODRecoDecayHF3Prong *part,
2680                                                             AliRDHFCutsLctopKpi *cuts,
2681                                                             AliAODEvent* aod,
2682                                                             TClonesArray *arrMC)
2683 {
2684
2685   // FillAPrioriConcentrations
2686   cuts->SetUsePID(kFALSE); //Annalisa
2687   Int_t isSelected3ProngByLc=cuts->IsSelected(part,AliRDHFCuts::kCandidate,aod);
2688
2689   TString fillthis="";
2690
2691   if(isSelected3ProngByLc>0 && fReadMC) {
2692
2693     for (Int_t ii=0; ii<3; ii++) {
2694       AliAODTrack *prongTest=(AliAODTrack*)part->GetDaughter(ii);
2695       if (!prongTest) continue;
2696       Int_t labprongTest = prongTest->GetLabel();
2697       if(labprongTest<0) continue;
2698       AliAODMCParticle *mcprongTest = (AliAODMCParticle*)arrMC->At(labprongTest);
2699
2700       switch (TMath::Abs(mcprongTest->GetPdgCode())) {
2701       case 11:
2702         fillthis="hElIn3Prong";
2703         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2704         break;
2705       case 13:
2706         fillthis="hMuIn3Prong";
2707         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2708         break;
2709       case 211:
2710         fillthis="hPiIn3Prong";
2711         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2712         break;
2713       case 321:
2714         fillthis="hKaIn3Prong";
2715         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2716         break;
2717       case 2212:
2718         fillthis="hPrIn3Prong";
2719         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2720         break;
2721       default:
2722         break;
2723       }
2724     }
2725   }
2726
2727   cuts->SetUsePID(kTRUE); //Annalisa
2728
2729 }
2730
2731 //-----------------------
2732 Bool_t AliAnalysisTaskSELambdac::Is3ProngFromPDG(AliAODRecoDecayHF3Prong *d,
2733                                                     TClonesArray *arrayMC,
2734                                                     Int_t pdgToBeCompared)
2735 {
2736   //
2737   // Returs kTRUE if at least one of tracks in the current 3prong
2738   // cames from particle with pdgCode=pdgToBeCompared
2739   //
2740
2741   Bool_t localFlag = kFALSE;
2742
2743   for (Int_t ii=0;ii<3;ii++) {
2744     AliAODTrack *daugh = (AliAODTrack*)d->GetDaughter(ii);
2745
2746     Int_t lab = daugh->GetLabel();
2747     if(lab<0) continue;
2748     AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab);
2749     //Int_t partPdg = part->GetPdgCode();
2750     //printf(" -------------------------------------- partLab=%d partPdg=%d ---\n", lab, partPdg);
2751
2752     Int_t motherLabel = part->GetMother();
2753     if(motherLabel<0) continue;
2754
2755     AliAODMCParticle *motherPart = 0;
2756     Int_t motherPdg = 0;
2757     while (!localFlag && motherLabel>=0) {
2758       motherPart = (AliAODMCParticle*)arrayMC->At(motherLabel);
2759       motherPdg = motherPart->GetPdgCode();
2760       //printf(" ------------- motherLab=%d motherPdg=%d ---\n", motherLabel, motherPdg);
2761       if (TMath::Abs(motherPdg)==pdgToBeCompared) {
2762         //printf("-------------------------- trovato quark c/cbar\n");
2763         localFlag = kTRUE;
2764       }
2765       else {
2766         motherLabel = motherPart->GetMother();
2767       }
2768     }
2769
2770   }
2771
2772   return localFlag;
2773
2774 }
2775
2776 //-----------------------
2777 Bool_t AliAnalysisTaskSELambdac::IsTrackFromPDG(const AliAODTrack *daugh,
2778                                                    TClonesArray *arrayMC,
2779                                                    Int_t pdgToBeCompared)
2780 {
2781   //
2782   // Returs kTRUE if at the tracks comes
2783   // from particle with pdgCode=pdgToBeCompared
2784   //
2785
2786   Bool_t localFlag = kFALSE;
2787
2788   Int_t lab = daugh->GetLabel();
2789   if(lab<0) return localFlag;
2790   AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab);
2791   //Int_t partPdg = part->GetPdgCode();
2792   //printf(" -------------------------------------- partLab=%d partPdg=%d ---\n", lab, partPdg);
2793
2794   Int_t motherLabel = part->GetMother();
2795   if(motherLabel<0) return localFlag;
2796
2797   AliAODMCParticle *motherPart = 0;
2798   Int_t motherPdg = 0;
2799   while (!localFlag && motherLabel>=0) {
2800     motherPart = (AliAODMCParticle*)arrayMC->At(motherLabel);
2801     motherPdg = motherPart->GetPdgCode();
2802     //printf(" ------------- motherLab=%d motherPdg=%d ---\n", motherLabel, motherPdg);
2803     if (TMath::Abs(motherPdg)==pdgToBeCompared) {
2804       //printf("-------------------------- trovato quark c/cbar\n");
2805       localFlag = kTRUE;
2806     }
2807     else {
2808       motherLabel = motherPart->GetMother();
2809     }
2810   }
2811
2812   return localFlag;
2813
2814 }
2815
2816
2817
2818
2819 //-----------------------
2820 Bool_t AliAnalysisTaskSELambdac::IsThereAGeneratedLc(TClonesArray *arrayMC)
2821 {
2822   //
2823   // Returs kTRUE if there is a lepton related to the LambdaC
2824   //
2825
2826   Bool_t localFlag = kFALSE;
2827
2828   AliAODMCParticle *searchLc;
2829
2830   for (Int_t iii=0; iii<arrayMC->GetEntries(); iii++) {
2831     searchLc = (AliAODMCParticle*)arrayMC->At(iii);
2832     Int_t searchLcpdg =  searchLc->GetPdgCode();
2833     if (TMath::Abs(searchLcpdg) == 4122){
2834       localFlag = kTRUE;
2835       break;
2836     }
2837   }
2838
2839
2840   return localFlag;
2841
2842 }
2843 //_________________________________
2844 Int_t AliAnalysisTaskSELambdac::NumberPrimaries(const AliAODEvent *aods)
2845 {
2846 ///////////////estimate primaries
2847   Int_t counter=0;
2848
2849
2850   TClonesArray *aodtracks=(TClonesArray *)aods->GetTracks();
2851
2852   // for(Int_t ji=0;ji<aods->GetNTracks();ji++)
2853   for(Int_t ji=0;ji<aodtracks->GetEntriesFast();ji++)
2854     {
2855       AliAODTrack*aodTrack=(AliAODTrack*)aodtracks->UncheckedAt(ji);
2856       if(aodTrack->IsPrimaryCandidate()) counter++;
2857     }
2858   return counter;
2859
2860 }
2861 //_________________________________
2862
2863 void AliAnalysisTaskSELambdac::MultiplicityStudies(AliAODRecoDecayHF3Prong *part,
2864                                                       AliRDHFCutsLctopKpi *cuts,
2865                                                       AliAODEvent* aod,
2866                                                       TClonesArray *arrMC,
2867                                                       Bool_t &flag1,Bool_t &flag2,Bool_t &flag3,
2868                                                       Bool_t &flag4, Bool_t &flag5, Bool_t &flag6)
2869 {
2870
2871
2872 // Multiplicity studies
2873
2874   TString fillthis="";
2875
2876   Int_t pdgDgLctopKpi[3]={2212,321,211};
2877   Int_t lab=-9999;
2878   if(fReadMC)
2879     lab=part->MatchToMC(4122,arrMC,3,pdgDgLctopKpi); //return MC particle label if the array corresponds to a Lc, -1 if not (cf. AliAODRecoDecay.cxx)
2880
2881   cuts->SetUsePID(kFALSE); //Annalisa
2882   Int_t isSelected3ProngByLc=cuts->IsSelected(part,AliRDHFCuts::kCandidate,aod);
2883
2884   if(isSelected3ProngByLc>0 && fReadMC) {
2885     flag1 = kTRUE;
2886     if (lab>0) {
2887       flag3 = kTRUE;
2888     }
2889
2890     Bool_t is3ProngFromJPsi = Is3ProngFromPDG(part,arrMC,443);
2891     if (is3ProngFromJPsi) flag6=is3ProngFromJPsi;
2892
2893     Bool_t is3ProngFromC = Is3ProngFromPDG(part,arrMC,4);
2894     if (is3ProngFromC) flag4=is3ProngFromC;
2895
2896     Bool_t is3ProngFromB = Is3ProngFromPDG(part,arrMC,5);
2897     if (is3ProngFromB) flag5=is3ProngFromB;
2898
2899     for (Int_t ii=0; ii<3; ii++) {
2900       AliAODTrack *prongTest=(AliAODTrack*)part->GetDaughter(ii);
2901       if (!prongTest) continue;
2902       Int_t labprongTest = prongTest->GetLabel();
2903       if(labprongTest<0) continue;
2904       AliAODMCParticle *mcprongTest = (AliAODMCParticle*)arrMC->At(labprongTest);
2905
2906       switch (TMath::Abs(mcprongTest->GetPdgCode())) {
2907       case 11:
2908         fillthis="hElIn3Prong";
2909         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2910         if (IsTrackFromPDG(prongTest,arrMC,443)) {
2911           fillthis="hElIn3Prong6";
2912           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2913         }
2914         break;
2915       case 13:
2916         fillthis="hMuIn3Prong";
2917         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2918         if (IsTrackFromPDG(prongTest,arrMC,443)) {
2919           fillthis="hMuIn3Prong6";
2920           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2921         }
2922         break;
2923       case 211:
2924         fillthis="hPiIn3Prong";
2925         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2926         if (IsTrackFromPDG(prongTest,arrMC,443)) {
2927           fillthis="hPiIn3Prong6";
2928           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2929         }
2930         break;
2931       case 321:
2932         fillthis="hKaIn3Prong";
2933         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2934         if (IsTrackFromPDG(prongTest,arrMC,443)) {
2935           fillthis="hKaIn3Prong6";
2936           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2937         }
2938         break;
2939       case 2212:
2940         fillthis="hPrIn3Prong";
2941         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2942         if (IsTrackFromPDG(prongTest,arrMC,443)) {
2943           fillthis="hPrIn3Prong6";
2944           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2945         }
2946         break;
2947       default:
2948         break;
2949       }
2950
2951       if (lab>0) {
2952         switch (TMath::Abs(mcprongTest->GetPdgCode())) {
2953         case 11:
2954           fillthis="hElIn3Prong1";
2955           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2956           break;
2957         case 13:
2958           fillthis="hMuIn3Prong1";
2959           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2960           break;
2961         case 211:
2962           fillthis="hPiIn3Prong1";
2963           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2964           break;
2965         case 321:
2966           fillthis="hKaIn3Prong1";
2967           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2968           break;
2969         case 2212:
2970           fillthis="hPrIn3Prong1";
2971           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2972           break;
2973         default:
2974           break;
2975         }
2976       } else {
2977         switch (TMath::Abs(mcprongTest->GetPdgCode())) {
2978         case 11:
2979           fillthis="hElIn3Prong2";
2980           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2981           if (IsTrackFromPDG(prongTest,arrMC,4)) {
2982             fillthis="hElIn3Prong3";
2983             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2984           }
2985           if (IsTrackFromPDG(prongTest,arrMC,5)) {
2986             fillthis="hElIn3Prong4";
2987             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2988           }
2989           if (!IsTrackFromPDG(prongTest,arrMC,4) && !IsTrackFromPDG(prongTest,arrMC,5)) {
2990             fillthis="hElIn3Prong5";
2991             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2992           }
2993           break;
2994         case 13:
2995           fillthis="hMuIn3Prong2";
2996           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2997           if (IsTrackFromPDG(prongTest,arrMC,4)) {
2998             fillthis="hMuIn3Prong3";
2999             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3000           }
3001           if (IsTrackFromPDG(prongTest,arrMC,5)) {
3002             fillthis="hMuIn3Prong4";
3003             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3004           }
3005           if (!IsTrackFromPDG(prongTest,arrMC,4) && !IsTrackFromPDG(prongTest,arrMC,5)) {
3006             fillthis="hMuIn3Prong5";
3007             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3008           }
3009           break;
3010         case 211:
3011           fillthis="hPiIn3Prong2";
3012           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3013           if (IsTrackFromPDG(prongTest,arrMC,4)) {
3014             fillthis="hPiIn3Prong3";
3015             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3016           }
3017           if (IsTrackFromPDG(prongTest,arrMC,5)) {
3018             fillthis="hPiIn3Prong4";
3019             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3020           }
3021           if (!IsTrackFromPDG(prongTest,arrMC,4) && !IsTrackFromPDG(prongTest,arrMC,5)) {
3022             fillthis="hPiIn3Prong5";
3023             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3024           }
3025           break;
3026         case 321:
3027           fillthis="hKaIn3Prong2";
3028           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3029           if (IsTrackFromPDG(prongTest,arrMC,4)) {
3030             fillthis="hKaIn3Prong3";
3031             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3032           }
3033           if (IsTrackFromPDG(prongTest,arrMC,5)) {
3034             fillthis="hKaIn3Prong4";
3035             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3036           }
3037           if (!IsTrackFromPDG(prongTest,arrMC,4) && !IsTrackFromPDG(prongTest,arrMC,5)) {
3038             fillthis="hKaIn3Prong5";
3039             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3040           }
3041           break;
3042         case 2212:
3043           fillthis="hPrIn3Prong2";
3044           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3045           if (IsTrackFromPDG(prongTest,arrMC,4)) {
3046             fillthis="hPrIn3Prong3";
3047             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3048           }
3049           if (IsTrackFromPDG(prongTest,arrMC,5)) {
3050             fillthis="hPrIn3Prong4";
3051             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3052           }
3053           if (!IsTrackFromPDG(prongTest,arrMC,4) && !IsTrackFromPDG(prongTest,arrMC,5)) {
3054             fillthis="hPrIn3Prong5";
3055             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3056           }
3057           break;
3058         default:
3059           break;
3060         }
3061       }
3062       /*
3063         if (is3ProngFromC) {
3064         switch (TMath::Abs(mcprongTest->GetPdgCode())) {
3065         case 11:
3066         fillthis="hElIn3Prong3";
3067         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3068         break;
3069         case 13:
3070         fillthis="hMuIn3Prong3";
3071         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3072         break;
3073         case 211:
3074         fillthis="hPiIn3Prong3";
3075         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3076         break;
3077         case 321:
3078         fillthis="hKaIn3Prong3";
3079         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3080         break;
3081         case 2212:
3082         fillthis="hPrIn3Prong3";
3083         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3084         break;
3085         default:
3086         break;
3087         }
3088         } else { // !is3ProngFromC
3089         if (is3ProngFromB) {
3090         switch (TMath::Abs(mcprongTest->GetPdgCode())) {
3091         case 11:
3092         fillthis="hElIn3Prong4";
3093         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3094         break;
3095         case 13:
3096         fillthis="hMuIn3Prong4";
3097         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3098         break;
3099         case 211:
3100         fillthis="hPiIn3Prong4";
3101         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3102         break;
3103         case 321:
3104         fillthis="hKaIn3Prong4";
3105         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3106         break;
3107         case 2212:
3108         fillthis="hPrIn3Prong4";
3109         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3110         break;
3111         default:
3112         break;
3113         }
3114         } else {//!is3ProngFromB
3115         switch (TMath::Abs(mcprongTest->GetPdgCode())) {
3116         case 11:
3117         fillthis="hElIn3Prong5";
3118         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3119         break;
3120         case 13:
3121         fillthis="hMuIn3Prong5";
3122         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3123         break;
3124         case 211:
3125         fillthis="hPiIn3Prong5";
3126         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3127         break;
3128         case 321:
3129         fillthis="hKaIn3Prong5";
3130         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3131         break;
3132         case 2212:
3133         fillthis="hPrIn3Prong5";
3134         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3135         break;
3136         default:
3137         break;
3138         }
3139         }
3140         }
3141       */
3142     }
3143   }
3144
3145   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
3146   Double_t invmasscut=0.05;
3147
3148   Double_t minvLcpKpi = part->InvMassLcpKpi();
3149   Double_t minvLcpiKp = part->InvMassLcpiKp();
3150
3151   cuts->SetUsePID(kTRUE); //Annalisa
3152   Int_t isSelected3ProngByLcPID=cuts->IsSelected(part,AliRDHFCuts::kCandidate,aod);
3153   if (isSelected3ProngByLcPID>0) {
3154     if (TMath::Abs(minvLcpKpi-mPDG)<invmasscut || TMath::Abs(minvLcpiKp-mPDG)<invmasscut) {
3155       flag2 = kTRUE;
3156     }
3157   }
3158
3159
3160 }