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$ */
  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(),100,fLowmasslimit,fUpmasslimit);
 359     fMassHist[index]->Sumw2();
 360     hisname.Form("hMassPt%dTC",i);
 361     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 362     fMassHistTC[index]->Sumw2();
 363
 364     hisname.Form("hMassPtLpi%d",i);
 365     fMassHistLpi[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 366     fMassHistLpi[index]->Sumw2();
 367     hisname.Form("hMassPtLpi%dTC",i);
 368     fMassHistLpiTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 369     fMassHistLpiTC[index]->Sumw2();
 370
 371     hisname.Form("hMassPtKp%d",i);
 372     fMassHistKp[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 373     fMassHistKp[index]->Sumw2();
 374     hisname.Form("hMassPtKp%dTC",i);
 375     fMassHistKpTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 376     fMassHistKpTC[index]->Sumw2();
 377     hisname.Form("hMassPtDk%d",i);
 378     fMassHistDk[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 379     fMassHistDk[index]->Sumw2();
 380     hisname.Form("hMassPtDk%dTC",i);
 381     fMassHistDkTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 382     fMassHistDkTC[index]->Sumw2();
 383
 384     hisname.Form("hMassPt3Pr%d",i);
 385     fMassHist3Pr[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 386     fMassHist3Pr[index]->Sumw2();
 387     hisname.Form("hMassPt3Pr%dTC",i);
 388     fMassHist3PrTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,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(),100,fLowmasslimit,fUpmasslimit);
 394     fMassHist[index]->Sumw2();
 395     hisname.Form("hSigPt%dTC",i);
 396     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 397     fMassHistTC[index]->Sumw2();
 398     hisname.Form("hSigPtLpi%d",i);
 399     fMassHistLpi[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 400     fMassHistLpi[index]->Sumw2();
 401     hisname.Form("hSigPtLpi%dTC",i);
 402     fMassHistLpiTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 403     fMassHistLpiTC[index]->Sumw2();
 404
 405     hisname.Form("hSigPtKp%d",i);
 406     fMassHistKp[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 407     fMassHistKp[index]->Sumw2();
 408     hisname.Form("hSigPtKp%dTC",i);
 409     fMassHistKpTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 410     fMassHistKpTC[index]->Sumw2();
 411
 412     hisname.Form("hSigPtDk%d",i);
 413     fMassHistDk[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 414     fMassHistDk[index]->Sumw2();
 415     hisname.Form("hSigPtDk%dTC",i);
 416     fMassHistDkTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 417     fMassHistDkTC[index]->Sumw2();
 418
 419     hisname.Form("hSigPt3Pr%d",i);
 420     fMassHist3Pr[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 421     fMassHist3Pr[index]->Sumw2();
 422     hisname.Form("hSigPt3Pr%dTC",i);
 423     fMassHist3PrTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,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(),100,fLowmasslimit,fUpmasslimit);
 429     fMassHist[index]->Sumw2();
 430     hisname.Form("hBkgPt%dTC",i);
 431     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 432     fMassHistTC[index]->Sumw2();
 433     hisname.Form("hBkgPtLpi%d",i);
 434     fMassHistLpi[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 435     fMassHistLpi[index]->Sumw2();
 436     hisname.Form("hBkgPtLpi%dTC",i);
 437     fMassHistLpiTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 438     fMassHistLpiTC[index]->Sumw2();
 439
 440     hisname.Form("hBkgPtKp%d",i);
 441     fMassHistKp[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 442     fMassHistKp[index]->Sumw2();
 443     hisname.Form("hBkgPtKp%dTC",i);
 444     fMassHistKpTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 445     fMassHistKpTC[index]->Sumw2();
 446
 447     hisname.Form("hBkgPtDk%d",i);
 448     fMassHistDk[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 449     fMassHistDk[index]->Sumw2();
 450     hisname.Form("hBkgPtDk%dTC",i);
 451     fMassHistDkTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 452     fMassHistDkTC[index]->Sumw2();
 453
 454     hisname.Form("hBkgPt3Pr%d",i);
 455     fMassHist3Pr[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 456     fMassHist3Pr[index]->Sumw2();
 457     hisname.Form("hBkgPt3Pr%dTC",i);
 458     fMassHist3PrTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,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",100,fLowmasslimit,fUpmasslimit);
 485   fhMassPtGreater3->Sumw2();
 486   fOutput->Add(fhMassPtGreater3);
 487   fhMassPtGreater3TC=new TH1F("fhMassPtGreater3TC","Pt > 3 GeV/c",100,fLowmasslimit,fUpmasslimit);
 488   fhMassPtGreater3TC->Sumw2();
 489   fOutput->Add(fhMassPtGreater3TC);
 490   fhMassPtGreater3Kp=new TH1F("fhMassPtGreater3Kp","Pt > 3 GeV/c",100,fLowmasslimit,fUpmasslimit);
 491   fhMassPtGreater3Kp->Sumw2();
 492   fOutput->Add(fhMassPtGreater3Kp);
 493   fhMassPtGreater3KpTC=new TH1F("fhMassPtGreater3KpTC","Pt > 3 GeV/c",100,fLowmasslimit,fUpmasslimit);
 494   fhMassPtGreater3KpTC->Sumw2();
 495   fOutput->Add(fhMassPtGreater3KpTC);
 496   fhMassPtGreater3Lpi=new TH1F("fhMassPtGreater3Lpi","Pt > 3 GeV/c",100,fLowmasslimit,fUpmasslimit);
 497   fhMassPtGreater3Lpi->Sumw2();
 498   fOutput->Add(fhMassPtGreater3Lpi);
 499   fhMassPtGreater3LpiTC=new TH1F("fhMassPtGreater3LpiTC","Pt > 3 GeV/c",100,fLowmasslimit,fUpmasslimit);
 500   fhMassPtGreater3LpiTC->Sumw2();
 501   fOutput->Add(fhMassPtGreater3LpiTC);
 502   fhMassPtGreater3Dk=new TH1F("fhMassPtGreater3Dk","Pt > 3 GeV/c",100,fLowmasslimit,fUpmasslimit);
 503   fhMassPtGreater3Dk->Sumw2();
 504   fOutput->Add(fhMassPtGreater3Dk);
 505   fhMassPtGreater3DkTC=new TH1F("fhMassPtGreater3DkTC","Pt > 3 GeV/c",100,fLowmasslimit,fUpmasslimit);
 506   fhMassPtGreater3DkTC->Sumw2();
 507   fOutput->Add(fhMassPtGreater3DkTC);
 508
 509   fhMassPtGreater33Pr=new TH1F("fhMassPtGreater33Pr","Pt > 3 GeV/c",100,fLowmasslimit,fUpmasslimit);
 510   fhMassPtGreater33Pr->Sumw2();
 511   fOutput->Add(fhMassPtGreater33Pr);
 512   fhMassPtGreater33PrTC=new TH1F("fhMassPtGreater33PrTC","Pt > 3 GeV/c",100,fLowmasslimit,fUpmasslimit);
 513   fhMassPtGreater33PrTC->Sumw2();
 514   fOutput->Add(fhMassPtGreater33PrTC);
 515
 516
 517   fhMassPtGreater2=new TH1F("fhMassPtGreater2","Pt > 2 GeV/c",100,fLowmasslimit,fUpmasslimit);
 518   fhMassPtGreater2->Sumw2();
 519   fOutput->Add(fhMassPtGreater2);
 520   fhMassPtGreater2TC=new TH1F("fhMassPtGreater2TC","Pt > 2 GeV/c",100,fLowmasslimit,fUpmasslimit);
 521   fhMassPtGreater2TC->Sumw2();
 522   fOutput->Add(fhMassPtGreater2TC);
 523   fhMassPtGreater2Kp=new TH1F("fhMassPtGreater2Kp","Pt > 2 GeV/c",100,fLowmasslimit,fUpmasslimit);
 524   fhMassPtGreater2Kp->Sumw2();
 525   fOutput->Add(fhMassPtGreater2Kp);
 526   fhMassPtGreater2KpTC=new TH1F("fhMassPtGreater2KpTC","Pt > 2 GeV/c",100,fLowmasslimit,fUpmasslimit);
 527   fhMassPtGreater2KpTC->Sumw2();
 528   fOutput->Add(fhMassPtGreater2KpTC);
 529   fhMassPtGreater2Lpi=new TH1F("fhMassPtGreater2Lpi","Pt > 2 GeV/c",100,fLowmasslimit,fUpmasslimit);
 530   fhMassPtGreater2Lpi->Sumw2();
 531   fOutput->Add(fhMassPtGreater2Lpi);
 532   fhMassPtGreater2LpiTC=new TH1F("fhMassPtGreater2LpiTC","Pt > 2 GeV/c",100,fLowmasslimit,fUpmasslimit);
 533   fhMassPtGreater2LpiTC->Sumw2();
 534   fOutput->Add(fhMassPtGreater2LpiTC);
 535   fhMassPtGreater2Dk=new TH1F("fhMassPtGreater2Dk","Pt > 2 GeV/c",100,fLowmasslimit,fUpmasslimit);
 536   fhMassPtGreater2Dk->Sumw2();
 537   fOutput->Add(fhMassPtGreater2Dk);
 538   fhMassPtGreater2DkTC=new TH1F("fhMassPtGreater2DkTC","Pt > 2 GeV/c",100,fLowmasslimit,fUpmasslimit);
 539   fhMassPtGreater2DkTC->Sumw2();
 540   fOutput->Add(fhMassPtGreater2DkTC);
 541   fhMassPtGreater23Pr=new TH1F("fhMassPtGreater23Pr","Pt > 2 GeV/c",100,fLowmasslimit,fUpmasslimit);
 542   fhMassPtGreater23Pr->Sumw2();
 543   fOutput->Add(fhMassPtGreater23Pr);
 544   fhMassPtGreater23PrTC=new TH1F("fhMassPtGreater23PrTC","Pt > 2 GeV/c",100,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(),100,fLowmasslimit,fUpmasslimit);
 572   fOutputMC->Add(hMassInv);
 573   hisname.Form("hbMass");
 574   TH1F *hBMassInv=new TH1F(hisname.Data(),hisname.Data(),100,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   fCounter = new AliNormalizationCounter("NormalizationCounter");
1072    fCounter->Init();
1073     PostData(7,fCounter);
1074   if (fFillNtuple) {
1075     //OpenFile(3); // 2 is the slot number of the ntuple
1076     fNtupleLambdac = new TNtuple("fNtupleLambdac","D +",
1077                                  "pdg:Px:Py:Pz:PtTrue:VxTrue:VyTrue:VzTrue:Ptpi:PtK:Ptpi2:PtRec:PointingAngle:DecLeng:VxRec:VyRec:VzRec:InvMass:sigvert:d0Pi:d0K:d0Pi2:dca:d0square");
1078     PostData(8,fNtupleLambdac);
1079   }
1080
1081   return;
1082 }
1083
1084 //________________________________________________________________________
1085 void AliAnalysisTaskSELambdac::UserExec(Option_t */*option*/)
1086 {
1087   // Execute analysis for current event:
1088   // heavy flavor candidates association to MC truth
1089
1090   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
1091   //tmp
1092   fHistNEvents->Fill(0); // count event
1093   // Post the data already here
1094
1095   TClonesArray *array3Prong = 0;
1096   TClonesArray *arrayLikeSign =0;
1097   if(!aod && AODEvent() && IsStandardAOD()) {
1098     // In case there is an AOD handler writing a standard AOD, use the AOD
1099     // event in memory rather than the input (ESD) event.
1100     aod = dynamic_cast<AliAODEvent*> (AODEvent());
1101     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
1102     // have to taken from the AOD event hold by the AliAODExtension
1103     AliAODHandler* aodHandler = (AliAODHandler*)
1104       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
1105     if(aodHandler->GetExtensions()) {
1106       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
1107       AliAODEvent *aodFromExt = ext->GetAOD();
1108       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
1109       arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
1110     }
1111   } else if(aod) {
1112     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
1113     arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
1114   }
1115
1116   if(!aod) return;
1117
1118   TClonesArray *arrayMC=0;
1119   AliAODMCHeader *mcHeader=0;
1120
1121   // load MC particles
1122   if(fReadMC){
1123
1124     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1125     if(!arrayMC) {
1126       printf("AliAnalysisTaskSELambdac::UserExec: MC particles branch not found!\n");
1127       return;
1128     }
1129
1130
1131     // load MC header
1132     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1133     if(!mcHeader) {
1134       printf("AliAnalysisTaskSELambdac::UserExec: MC header branch not found!\n");
1135       return;
1136     }
1137   }
1138
1139   TString fillthis="";
1140   Int_t numberOfPrimaries= NumberPrimaries(aod);
1141
1142   if (fMultiplicityHists && fReadMC) {
1143     fillthis="hPrimariesvsAOD";
1144     ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(aod->GetNTracks(),numberOfPrimaries);
1145
1146     fillthis="hAll2MultiplicityInEvent";
1147     ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(aod->GetNTracks());
1148
1149     fillthis="hAll2MultiplicityPrimaryInEvent";
1150     ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1151
1152     if (IsThereAGeneratedLc(arrayMC)) {
1153       fillthis="h2MultiplicityInLcEvent";
1154       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1155     }
1156
1157   }
1158
1159   if(!array3Prong || !aod) {
1160     printf("AliAnalysisTaskSELambdac::UserExec: Charm3Prong branch not found!\n");
1161     return;
1162   }
1163   if(!arrayLikeSign) {
1164     printf("AliAnalysisTaskSELambdac::UserExec: LikeSign3Prong branch not found!\n");
1165     //  return;
1166   }
1167
1168   fNentries->Fill(0);
1169   fCounter->StoreEvent(aod,fRDCutsProduction,fReadMC);
1170   TString trigclass=aod->GetFiredTriggerClasses();
1171   if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fNentries->Fill(14);
1172   if(!fRDCutsProduction->IsEventSelected(aod)) {
1173     if(fRDCutsProduction->GetWhyRejection()==1) // rejected for pileup
1174       fNentries->Fill(13);
1175     return;
1176   }
1177
1178   // AOD primary vertex
1179   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
1180   if(!vtx1) return;
1181
1182   Bool_t isThereA3prongWithGoodTracks = kFALSE;
1183   Bool_t isThereA3ProngLcKine = kFALSE;
1184   Bool_t isThereA3ProngLcKineANDpid = kFALSE;
1185   Bool_t isThereA3ProngLcMC = kFALSE;
1186   Bool_t isThereA3ProngCyes = kFALSE;
1187   Bool_t isThereA3ProngByes = kFALSE;
1188   Bool_t isThereA3ProngJPsi = kFALSE;
1189
1190   Int_t n3Prong = array3Prong->GetEntriesFast();
1191   Int_t nSelectedloose[1]={0};
1192   Int_t nSelectedtight[1]={0};
1193
1194   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
1195     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
1196
1197     Bool_t unsetvtx=kFALSE;
1198     if(!d->GetOwnPrimaryVtx()){
1199       d->SetOwnPrimaryVtx(vtx1);
1200       unsetvtx=kTRUE;
1201     }
1202
1203     //if(d->GetSelectionMap()) {if(!d->HasSelectionBit(AliRDHFCuts::kLcCuts)) continue;}
1204     //if(d->GetSelectionMap()) {if(!d->HasSelectionBit(AliRDHFCuts::kLcPID)) continue;}
1205
1206     Int_t isSelectedTracks = fRDCutsProduction->IsSelected(d,AliRDHFCuts::kTracks,aod);
1207     if(!isSelectedTracks) continue;
1208
1209     isThereA3prongWithGoodTracks=kTRUE;
1210
1211     if (fRDCutsProduction->IsInFiducialAcceptance(d->Pt(),d->Y(4122))) fNentries->Fill(6);
1212
1213     Int_t ptbin=fRDCutsProduction->PtBin(d->Pt());
1214     if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
1215
1216     FillMassHists(aod,d,arrayMC,fRDCutsProduction,nSelectedloose,nSelectedtight);
1217     if(fFillVarHists) FillVarHists(d,arrayMC,fRDCutsProduction,/*fOutputMC,*/aod);
1218
1219     if (fPriorsHists && fReadMC)
1220       FillAPrioriConcentrations(d, fRDCutsProduction, aod, arrayMC);
1221     if (fMultiplicityHists && fReadMC)
1222       MultiplicityStudies(d, fRDCutsProduction, aod, arrayMC,
1223                           isThereA3ProngLcKine,isThereA3ProngLcKineANDpid,isThereA3ProngLcMC,
1224                           isThereA3ProngCyes,isThereA3ProngByes,isThereA3ProngJPsi);
1225
1226     /*
1227     //start OS analysis
1228     if(labDp<0)fHistOSbkg->Fill(d->InvMassDplus());
1229     fHistOS->Fill(d->InvMassDplus());
1230     */
1231
1232     if(unsetvtx) d->UnsetOwnPrimaryVtx();
1233   }
1234   fCounter->StoreCandidates(aod,nSelectedloose[0],kTRUE);
1235   fCounter->StoreCandidates(aod,nSelectedtight[0],kFALSE);
1236
1237   if (fMultiplicityHists && fReadMC) {
1238     fillthis="hAllMultiplicityInEvent";
1239     ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(aod->GetNTracks());
1240
1241     fillthis="hAllMultiplicityPrimaryInEvent";
1242     ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1243
1244     if (IsThereAGeneratedLc(arrayMC)) {
1245       fillthis="hMultiplicityInLcEvent";
1246       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1247     }
1248
1249     if (isThereA3prongWithGoodTracks) {
1250       fillthis="hMultiplicityInEvent";
1251       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1252     }
1253     if (isThereA3ProngLcKine) {
1254       fillthis="hMultiplicityIn3ProngLC";
1255       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1256     }
1257     if (isThereA3ProngLcKineANDpid) {
1258       fillthis="hMultiplicityInLCpid";
1259       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1260     }
1261     if (isThereA3ProngLcMC) {
1262       fillthis="hMultiplicityInLCmc";
1263       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1264     }
1265     if (isThereA3ProngLcKine && !isThereA3ProngLcMC) {
1266       fillthis="hMultiplicityInLCNomc";
1267       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1268     }
1269
1270     if (isThereA3ProngCyes) {
1271       fillthis="hMultiplicityYesC";
1272       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1273     }
1274     if (isThereA3ProngByes) {
1275       fillthis="hMultiplicityYesB";
1276       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1277     }
1278     if (isThereA3ProngJPsi) {
1279       fillthis="hMultiplicityJPsi";
1280       ((TH1I*)fMultiplicity->FindObject(fillthis))->Fill(numberOfPrimaries);
1281     }
1282
1283   }
1284
1285
1286   PostData(1,fOutput);
1287   if (fFillVarHists) PostData(3,fOutputMC);
1288   if (fPriorsHists) PostData(5,fAPriori);
1289   if (fMultiplicityHists) PostData(6,fMultiplicity);
1290
1291   PostData(4,fNentries);
1292   PostData(7,fCounter);
1293
1294   return;
1295 }
1296
1297
1298
1299 //________________________________________________________________________
1300 void AliAnalysisTaskSELambdac::Terminate(Option_t */*option*/)
1301 {
1302   // Terminate analysis
1303   //
1304   if (fDebug > 1) printf("AnalysisTaskSELambdac: Terminate() \n");
1305
1306   fOutput = dynamic_cast<TList*> (GetOutputData(1));
1307   if (!fOutput) {
1308     printf("ERROR: fOutput not available\n");
1309     return;
1310   }
1311   //fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
1312
1313   if(fFillNtuple){
1314     fNtupleLambdac = dynamic_cast<TNtuple*>(GetOutputData(3));
1315   }
1316
1317
1318   return;
1319 }
1320
1321 //________________________________________________________________________
1322 Int_t AliAnalysisTaskSELambdac::MatchToMCLambdac(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC) const{
1323   // check if the candidate is a Lambdac decaying in pKpi or in the resonant channels
1324   Int_t lambdacLab[3]={0,0,0};
1325   Int_t pdgs[3]={0,0,0};
1326   for(Int_t i=0;i<3;i++){
1327     AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
1328     Int_t lab=daugh->GetLabel();
1329     if(lab<0) return 0;
1330     AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab);
1331     if(!part) continue;
1332     pdgs[i]=part->GetPdgCode();
1333     Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
1334     if(partPdgcode==211 || partPdgcode==321 || partPdgcode==2212){
1335       Int_t motherLabel=part->GetMother();
1336       if(motherLabel<0) return 0;
1337       AliAODMCParticle *motherPart = (AliAODMCParticle*)arrayMC->At(motherLabel);
1338       if(!motherPart) continue;
1339       Int_t motherPdg = TMath::Abs(motherPart->GetPdgCode());
1340       if(motherPdg==4122) {
1341         if(GetLambdacDaugh(motherPart,arrayMC)){lambdacLab[i]=motherLabel;continue;}
1342       }
1343       if(motherPdg==313 || motherPdg==2224 || motherPdg==3124){
1344         Int_t granMotherLabel=motherPart->GetMother();
1345         if(granMotherLabel<0) return 0;
1346         AliAODMCParticle *granMotherPart = (AliAODMCParticle*)arrayMC->At(granMotherLabel);
1347         if(!granMotherPart) continue;
1348         Int_t granMotherPdg  = TMath::Abs(granMotherPart->GetPdgCode());
1349         if(granMotherPdg ==4122) {
1350           if(GetLambdacDaugh(granMotherPart,arrayMC)) {lambdacLab[i]=granMotherLabel;continue;}
1351         }
1352       }
1353     }
1354   }
1355
1356   if(lambdacLab[0]==lambdacLab[1] && lambdacLab[1]==lambdacLab[2]) {return lambdacLab[0];}
1357   return 0;
1358
1359 }
1360 //------------------------
1361 Bool_t AliAnalysisTaskSELambdac::GetLambdacDaugh(AliAODMCParticle *part,TClonesArray *arrayMC) const{
1362   // check if the particle is a lambdac and if its decay mode is the correct one
1363   Int_t numberOfLambdac=0;
1364   if(TMath::Abs(part->GetPdgCode())!=4122) return kFALSE;
1365   Int_t daughTmp[2];
1366   daughTmp[0]=part->GetDaughter(0);
1367   daughTmp[1]=part->GetDaughter(1);
1368   Int_t nDaugh = (Int_t)part->GetNDaughters();
1369   if(nDaugh<2) return kFALSE;
1370   if(nDaugh>3) return kFALSE;
1371   AliAODMCParticle* pdaugh1 = (AliAODMCParticle*)arrayMC->At(part->GetDaughter(0));
1372   if(!pdaugh1) {return kFALSE;}
1373   Int_t number1 = TMath::Abs(pdaugh1->GetPdgCode());
1374   AliAODMCParticle* pdaugh2 = (AliAODMCParticle*)arrayMC->At(part->GetDaughter(1));
1375   if(!pdaugh2) {return kFALSE;}
1376   Int_t number2 = TMath::Abs(pdaugh2->GetPdgCode());
1377
1378   if(nDaugh==3){
1379     Int_t thirdDaugh=part->GetDaughter(1)-1;
1380     AliAODMCParticle* pdaugh3 = (AliAODMCParticle*)arrayMC->At(thirdDaugh);
1381     Int_t number3 = TMath::Abs(pdaugh3->GetPdgCode());
1382     if((number1==321 && number2==211 && number3==2212) ||
1383        (number1==211 && number2==321 && number3==2212) ||
1384        (number1==211 && number2==2212 && number3==321) ||
1385        (number1==321 && number2==2212 && number3==211) ||
1386        (number1==2212 && number2==321 && number3==211) ||
1387        (number1==2212 && number2==211 && number3==321)) {
1388       numberOfLambdac++;
1389     }
1390   }
1391
1392   if(nDaugh==2){
1393
1394     //Lambda resonant
1395
1396     //Lambda -> p K*0
1397     //
1398     Int_t nfiglieK=0;
1399
1400     if((number1==2212 && number2==313)){
1401       nfiglieK=pdaugh2->GetNDaughters();
1402       if(nfiglieK!=2) return kFALSE;
1403       AliAODMCParticle* pdaughK1 = (AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
1404       AliAODMCParticle* pdaughK2 = (AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
1405       if(!pdaughK1) return kFALSE;
1406       if(!pdaughK2) return kFALSE;
1407       if((TMath::Abs(pdaughK1->GetPdgCode())==211 && TMath::Abs(pdaughK2->GetPdgCode())==321) || (TMath::Abs(pdaughK1->GetPdgCode())==321 && TMath::Abs(pdaughK2->GetPdgCode())==211)) numberOfLambdac++;
1408     }
1409
1410     if((number1==313 && number2==2212)){
1411       nfiglieK=pdaugh1->GetNDaughters();
1412       if(nfiglieK!=2) return kFALSE;
1413       AliAODMCParticle* pdaughK1 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
1414       AliAODMCParticle* pdaughK2 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
1415       if(!pdaughK1) return kFALSE;
1416       if(!pdaughK2) return kFALSE;
1417       if((TMath::Abs(pdaughK1->GetPdgCode())==211 && TMath::Abs(pdaughK2->GetPdgCode())==321) || (TMath::Abs(pdaughK1->GetPdgCode())==321 && TMath::Abs(pdaughK2->GetPdgCode())==211)) numberOfLambdac++;
1418     }
1419
1420     //Lambda -> Delta++ k
1421     Int_t nfiglieDelta=0;
1422     if(number1==321 && number2==2224){
1423       nfiglieDelta=pdaugh2->GetNDaughters();
1424       if(nfiglieDelta!=2) return kFALSE;
1425       AliAODMCParticle *pdaughD1=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
1426       AliAODMCParticle *pdaughD2=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
1427       if(!pdaughD1) return kFALSE;
1428       if(!pdaughD2) return kFALSE;
1429       if((TMath::Abs(pdaughD1->GetPdgCode())==211 && TMath::Abs(pdaughD2->GetPdgCode())==2212) || (TMath::Abs(pdaughD1->GetPdgCode())==2212 && TMath::Abs(pdaughD2->GetPdgCode())==211)) numberOfLambdac++;
1430     }
1431     if(number1==2224 && number2==321){
1432       nfiglieDelta=pdaugh1->GetNDaughters();
1433       if(nfiglieDelta!=2) return kFALSE;
1434       AliAODMCParticle* pdaughD1 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
1435       AliAODMCParticle* pdaughD2 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
1436       if(!pdaughD1) return kFALSE;
1437       if(!pdaughD2) return kFALSE;
1438       if((TMath::Abs(pdaughD1->GetPdgCode())==211 && TMath::Abs(pdaughD2->GetPdgCode())==2212) || (TMath::Abs(pdaughD1->GetPdgCode())==2212 && TMath::Abs(pdaughD2->GetPdgCode())==211)) numberOfLambdac++;
1439     }
1440
1441
1442     //Lambdac -> Lambda(1520) pi
1443     Int_t nfiglieLa=0;
1444     if(number1==3124 && number2==211){
1445       nfiglieLa=pdaugh1->GetNDaughters();
1446       if(nfiglieLa!=2) return kFALSE;
1447       AliAODMCParticle *pdaughL1=(AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
1448       AliAODMCParticle *pdaughL2=(AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
1449       if(!pdaughL1) return kFALSE;
1450       if(!pdaughL2) return kFALSE;
1451       if((TMath::Abs(pdaughL1->GetPdgCode())==321 && TMath::Abs(pdaughL2->GetPdgCode())==2212) || (TMath::Abs(pdaughL1->GetPdgCode())==2212 && TMath::Abs(pdaughL2->GetPdgCode())==321)) numberOfLambdac++;
1452     }
1453     if(number1==211 && number2==3124){
1454       nfiglieLa=pdaugh2->GetNDaughters();
1455       if(nfiglieLa!=2) return kFALSE;
1456       AliAODMCParticle *pdaughL1=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
1457       AliAODMCParticle *pdaughL2=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
1458       if(!pdaughL1) return kFALSE;
1459       if(!pdaughL2) return kFALSE;
1460       if((TMath::Abs(pdaughL1->GetPdgCode())==321 && TMath::Abs(pdaughL2->GetPdgCode())==2212) || (TMath::Abs(pdaughL1->GetPdgCode())==2212 && TMath::Abs(pdaughL2->GetPdgCode())==321)) numberOfLambdac++;
1461
1462     }
1463   }
1464
1465   if(numberOfLambdac>0) {return kTRUE;}
1466   return kFALSE;
1467 }
1468 //-----------------------------
1469 Bool_t AliAnalysisTaskSELambdac::IspKpiMC(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC) const{
1470   // Apply MC PID
1471   Int_t lab[3]={0,0,0},pdgs[3]={0,0,0};
1472   for(Int_t i=0;i<3;i++){
1473     AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
1474     lab[i]=daugh->GetLabel();
1475     if(lab[i]<0) return kFALSE;
1476     AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab[i]);
1477     if(!part) return kFALSE;
1478     pdgs[i]=TMath::Abs(part->GetPdgCode());
1479   }
1480
1481   if(pdgs[0]==2212 && pdgs[1]==321 && pdgs[2]==211) return kTRUE;
1482
1483   return kFALSE;
1484 }
1485 //-----------------------------
1486 Bool_t AliAnalysisTaskSELambdac::IspiKpMC(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC) const{
1487
1488   // Apply MC PID
1489   Int_t lab[3]={0,0,0},pdgs[3]={0,0,0};
1490   for(Int_t i=0;i<3;i++){
1491     AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
1492     lab[i]=daugh->GetLabel();
1493     if(lab[i]<0) return kFALSE;
1494     AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab[i]);
1495     if(!part) return kFALSE;
1496     pdgs[i]=TMath::Abs(part->GetPdgCode());
1497   }
1498
1499   if(pdgs[2]==2212 && pdgs[1]==321 && pdgs[0]==211) {return kTRUE;}
1500
1501   return kFALSE;
1502 }
1503 //--------------------------------------
1504 Bool_t AliAnalysisTaskSELambdac::VertexingKF(AliAODRecoDecayHF3Prong *d,Int_t *pdgs,Double_t field) const{
1505   // apply vertexing KF
1506   Int_t iprongs[3]={0,1,2};
1507   Double_t mass[2]={0.,0.};
1508   Bool_t constraint=kFALSE;
1509   if(fCutsKF[0]>0.)constraint=kTRUE;
1510   //topological constr
1511   AliKFParticle *lambdac=d->ApplyVertexingKF(iprongs,3,pdgs,constraint,field,mass);
1512   if(!lambdac) return kFALSE;
1513   if(lambdac->GetChi2()/lambdac->GetNDF()>fCutsKF[1]) return kFALSE;
1514   return kTRUE;
1515 }
1516 //-------------------------------------
1517 void AliAnalysisTaskSELambdac::IspiKpResonant(AliAODRecoDecayHF3Prong *d,Double_t field,Int_t *resNumber) const{
1518
1519   // apply PID using the resonant channels
1520   //if lambda* -> pk
1521   Double_t mass[2]={1.520,0.005};
1522   Int_t ipRes[2]={1,2};
1523   Int_t pdgres[2]={321,2212};
1524   AliKFParticle *lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1525   Double_t probLa=TMath::Prob(lambda1520->GetChi2(),lambda1520->GetNDF());
1526   if(probLa>0.1) resNumber[0]=1;
1527   //K* -> kpi
1528   mass[0]=0.8961;mass[1]=0.03;
1529   ipRes[0]=0;ipRes[1]=1;
1530   pdgres[0]=211;pdgres[1]=321;
1531   AliKFParticle *kstar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1532   Double_t probKa=TMath::Prob(kstar->GetChi2(),kstar->GetNDF());
1533   if(probKa>0.1) resNumber[1]=1;
1534
1535  // Delta++
1536   mass[0]=1.232;mass[1]=0.15;
1537   ipRes[0]=0;ipRes[1]=2;
1538   pdgres[0]=211;pdgres[1]=2122;
1539   AliKFParticle *delta=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1540   Double_t probDe=TMath::Prob(delta->GetChi2(),delta->GetNDF());
1541   if(probDe>0.1) resNumber[2]=1;
1542
1543   return ;
1544
1545 }
1546 //-------------------------------------
1547 void AliAnalysisTaskSELambdac::IspKpiResonant(AliAODRecoDecayHF3Prong *d,Double_t field,Int_t *resNumber) const{
1548
1549   // apply PID using the resonant channels
1550   //if lambda* -> pk
1551   Double_t mass[2]={1.520,0.005};
1552   Int_t ipRes[2]={0,1};
1553   Int_t pdgres[2]={2212,321};
1554   AliKFParticle *lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1555   Double_t probLa=TMath::Prob(lambda1520->GetChi2(),lambda1520->GetNDF());
1556   if(probLa>0.1) resNumber[0]=1;
1557   //K* -> kpi
1558   mass[0]=0.8961;mass[1]=0.03;
1559   ipRes[0]=1;ipRes[1]=2;
1560   pdgres[1]=211;pdgres[0]=321;
1561   AliKFParticle *kstar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1562   Double_t probKa=TMath::Prob(kstar->GetChi2(),kstar->GetNDF());
1563   if(probKa>0.1) resNumber[1]=1;
1564
1565   //    Delta++
1566   mass[0]=1.232;mass[1]=0.15;
1567   ipRes[0]=0;ipRes[1]=2;
1568   pdgres[0]=2122;pdgres[1]=211;
1569   AliKFParticle *delta=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1570   Double_t probDe=TMath::Prob(delta->GetChi2(),delta->GetNDF());
1571   if(probDe>0.1) resNumber[2]=1;
1572
1573   return ;
1574
1575 }
1576 //---------------------------
1577 void AliAnalysisTaskSELambdac::FillMassHists(AliAODEvent *aod,AliAODRecoDecayHF3Prong *part,
1578                                                 TClonesArray *arrayMC, AliRDHFCutsLctopKpi *cuts,Int_t *nSelectedloose,Int_t *nSelectedtight)
1579 {
1580
1581   //if MC PID or no PID, unset PID
1582   if(!fRealPid) cuts->SetUsePID(kFALSE);
1583   Int_t selection=cuts->IsSelected(part,AliRDHFCuts::kCandidate,aod);
1584   if(selection>0){
1585     nSelectedloose[0]=nSelectedloose[0]+1;
1586     Int_t iPtBin = -1;
1587     Double_t ptCand = part->Pt();
1588     Int_t index=0;
1589
1590     for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
1591       if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
1592     }
1593
1594     Float_t pdgCode=-2;
1595     Float_t pdgCode1=-1;
1596     Float_t pdgCode2=-1;
1597     Int_t labDp=-1;
1598     Float_t deltaPx=0.;
1599     Float_t deltaPy=0.;
1600     Float_t deltaPz=0.;
1601     Float_t truePt=0.;
1602     Float_t xDecay=0.;
1603     Float_t yDecay=0.;
1604     Float_t zDecay=0.;
1605
1606     if(fReadMC){
1607       labDp = MatchToMCLambdac(part,arrayMC);
1608       if(labDp>0){
1609         AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
1610         AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
1611         AliAODMCParticle *dg1 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(1));
1612         deltaPx=partDp->Px()-part->Px();
1613         deltaPy=partDp->Py()-part->Py();
1614         deltaPz=partDp->Pz()-part->Pz();
1615         truePt=partDp->Pt();
1616         xDecay=dg0->Xv();
1617         yDecay=dg0->Yv();
1618         zDecay=dg0->Zv();
1619         pdgCode=TMath::Abs(partDp->GetPdgCode());
1620         pdgCode1=TMath::Abs(dg0->GetPdgCode());
1621         pdgCode2=TMath::Abs(dg1->GetPdgCode());
1622
1623       }else{
1624         pdgCode=-1;
1625       }
1626     }
1627
1628     Double_t invMasspKpi=-1.;
1629     Double_t invMasspiKp=-1.;
1630     Double_t invMasspKpiLpi=-1.;
1631     Double_t invMasspiKpLpi=-1.;
1632     Double_t invMasspKpiKp=-1.;
1633     Double_t invMasspiKpKp=-1.;
1634     Double_t invMasspiKpDk=-1.;
1635     Double_t invMasspKpiDk=-1.;
1636     Double_t invMasspiKp3Pr=-1.;
1637     Double_t invMasspKpi3Pr=-1.;
1638     Double_t field=aod->GetMagneticField();
1639     //apply MC PID
1640     if(fReadMC && fMCPid){
1641       if(IspKpiMC(part,arrayMC)) invMasspKpi=part->InvMassLcpKpi();
1642       if(IspiKpMC(part,arrayMC)) invMasspiKp=part->InvMassLcpiKp();
1643     }
1644
1645     // apply realistic PID
1646     if(fRealPid){
1647       if(selection==1 || selection==3) invMasspKpi=part->InvMassLcpKpi();
1648       if(selection>=2) invMasspiKp=part->InvMassLcpiKp();
1649     }
1650     //apply PID using resonances
1651     if (fResPid && fRealPid) {
1652       Int_t ispKpi[3]={0,0,0};
1653       IspKpiResonant(part,field,ispKpi);
1654       Int_t ispiKp[3]={0,0,0};
1655       IspiKpResonant(part,field,ispiKp);
1656       if (selection==3 || selection==1) {
1657        if(ispKpi[0]==1){
1658         invMasspKpiLpi=part->InvMassLcpKpi();
1659        }
1660        if(ispKpi[1]==1){
1661         invMasspKpiKp=part->InvMassLcpKpi();
1662        }
1663        if(ispKpi[2]==1){
1664         invMasspKpiDk=part->InvMassLcpKpi();
1665        }
1666        if(ispKpi[2]==0 && ispKpi[1]==0 && ispKpi[0]==0){
1667         invMasspKpi3Pr=part->InvMassLcpKpi();
1668        }
1669       }
1670       if(selection>=2) {
1671        if(ispiKp[0]==1){
1672         invMasspiKpLpi=part->InvMassLcpiKp();
1673        }
1674        if(ispiKp[1]==1){
1675         invMasspiKpKp=part->InvMassLcpiKp();
1676        }
1677        if(ispiKp[2]==1){
1678         invMasspiKpDk=part->InvMassLcpiKp();
1679        }
1680        if(ispiKp[2]==0 && ispiKp[1]==0 && ispiKp[0]==0){
1681         invMasspiKp3Pr=part->InvMassLcpiKp();
1682        }
1683       }
1684
1685      }
1686
1687
1688     if(invMasspiKp<0. && invMasspKpi<0.) return;
1689
1690     Int_t passTightCuts=0;
1691     if(fAnalysis) {
1692      passTightCuts=fRDCutsAnalysis->IsSelected(part,AliRDHFCuts::kCandidate,aod);
1693      if(fUseKF){
1694       Int_t pdgs[3]={0,0,0};
1695       if(invMasspKpi>0.){
1696        pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
1697        if(!VertexingKF(part,pdgs,field)) {
1698         invMasspKpi=-1.;
1699         invMasspKpi3Pr=-1.;
1700         invMasspKpiDk=-1.;
1701         invMasspKpiLpi=-1.;
1702         invMasspKpiKp=-1.;
1703        }
1704       }
1705       if(invMasspiKp>0.){
1706         pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
1707         if(!VertexingKF(part,pdgs,field)) {
1708           invMasspiKp=-1.;
1709           invMasspiKp3Pr=-1.;
1710           invMasspiKpDk=-1.;
1711           invMasspiKpLpi=-1.;
1712           invMasspiKpKp=-1.;
1713        }
1714      }
1715    }
1716     if(passTightCuts>0) nSelectedtight[0]=nSelectedtight[0]+1;
1717   }
1718
1719     Float_t tmp[24];
1720     if (fFillNtuple) {
1721       tmp[0]=pdgCode;
1722       tmp[1]=deltaPx;
1723       tmp[2]=deltaPy;
1724       tmp[3]=deltaPz;
1725       tmp[4]=truePt;
1726       tmp[5]=xDecay;
1727       tmp[6]=yDecay;
1728       tmp[7]=zDecay;
1729       if(pdgCode1==2212) {tmp[8]=part->PtProng(0);}else{tmp[8]=0.;}
1730       if(pdgCode1==211) {tmp[10]=part->PtProng(0);}else{tmp[10]=0.;}
1731       tmp[9]=part->PtProng(1);
1732       if(pdgCode2==211) {tmp[10]=part->PtProng(2);}else{tmp[10]=0.;}
1733       tmp[11]=part->Pt();
1734       tmp[12]=part->CosPointingAngle();
1735       tmp[13]=part->DecayLength();
1736       tmp[14]=part->Xv();
1737       tmp[15]=part->Yv();
1738       tmp[16]=part->Zv();
1739       if(invMasspiKp>0.) tmp[17]=invMasspiKp;
1740       if(invMasspKpi>0.) tmp[17]=invMasspKpi;
1741       tmp[18]=part->GetSigmaVert();
1742       tmp[19]=part->Getd0Prong(0);
1743       tmp[20]=part->Getd0Prong(1);
1744       tmp[21]=part->Getd0Prong(2);
1745       tmp[22]=part->GetDCA();
1746       tmp[23]=part->Prodd0d0();
1747       fNtupleLambdac->Fill(tmp);
1748       PostData(7,fNtupleLambdac);
1749     }
1750
1751     if(part->Pt()>3.&& part->Pt()<=6.){
1752       if(invMasspiKp>0. && invMasspKpi>0.){
1753         if(invMasspiKp>0.) fhMassPtGreater3->Fill(invMasspiKp,0.5);
1754         if(invMasspKpi>0.) fhMassPtGreater3->Fill(invMasspKpi,0.5);
1755       }else{
1756         if(invMasspiKp>0.) fhMassPtGreater3->Fill(invMasspiKp);
1757         if(invMasspKpi>0.) fhMassPtGreater3->Fill(invMasspKpi);
1758       }
1759       if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1760         if(invMasspiKpLpi>0.) fhMassPtGreater3Lpi->Fill(invMasspiKpLpi,0.5);
1761         if(invMasspKpiLpi>0.) fhMassPtGreater3Lpi->Fill(invMasspKpiLpi,0.5);
1762       }else{
1763         if(invMasspiKpLpi>0.) fhMassPtGreater3Lpi->Fill(invMasspiKpLpi);
1764         if(invMasspKpiLpi>0.) fhMassPtGreater3Lpi->Fill(invMasspKpiLpi);
1765       }
1766       if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1767         if(invMasspiKpKp>0.) fhMassPtGreater3Kp->Fill(invMasspiKpKp,0.5);
1768         if(invMasspKpiKp>0.) fhMassPtGreater3Kp->Fill(invMasspKpiKp,0.5);
1769       }else{
1770         if(invMasspiKpKp>0.) fhMassPtGreater3Kp->Fill(invMasspiKpKp);
1771         if(invMasspKpiKp>0.) fhMassPtGreater3Kp->Fill(invMasspKpiKp);
1772       }
1773       if(invMasspiKpDk>0. && invMasspKpiDk>0.){
1774        if(invMasspiKpDk>0.) fhMassPtGreater3Dk->Fill(invMasspiKpDk,0.5);
1775        if(invMasspKpiDk>0.) fhMassPtGreater3Dk->Fill(invMasspKpiDk,0.5);
1776       }else{
1777        if(invMasspiKpDk>0.) fhMassPtGreater3Dk->Fill(invMasspiKpDk);
1778        if(invMasspKpiDk>0.) fhMassPtGreater3Dk->Fill(invMasspKpiDk);
1779       }
1780       if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0.){
1781         if(invMasspiKp3Pr>0.) fhMassPtGreater33Pr->Fill(invMasspiKp3Pr,0.5);
1782         if(invMasspKpi3Pr>0.) fhMassPtGreater33Pr->Fill(invMasspKpi3Pr,0.5);
1783       }else{
1784         if(invMasspiKp3Pr>0.) fhMassPtGreater33Pr->Fill(invMasspiKp3Pr);
1785         if(invMasspKpi3Pr>0.) fhMassPtGreater33Pr->Fill(invMasspKpi3Pr);
1786       }
1787
1788       if(passTightCuts>0){
1789         if(invMasspiKp>0. && invMasspKpi>0.){
1790           if(invMasspiKp>0.) fhMassPtGreater3TC->Fill(invMasspiKp,0.5);
1791           if(invMasspKpi>0.) fhMassPtGreater3TC->Fill(invMasspKpi,0.5);
1792         }else{
1793           if(invMasspiKp>0.) fhMassPtGreater3TC->Fill(invMasspiKp);
1794           if(invMasspKpi>0.) fhMassPtGreater3TC->Fill(invMasspKpi);
1795         }
1796         if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1797           if(invMasspiKpLpi>0.) fhMassPtGreater3LpiTC->Fill(invMasspiKpLpi,0.5);
1798           if(invMasspKpiLpi>0.) fhMassPtGreater3LpiTC->Fill(invMasspKpiLpi,0.5);
1799         }else{
1800           if(invMasspiKpLpi>0.) fhMassPtGreater3LpiTC->Fill(invMasspiKpLpi);
1801           if(invMasspKpiLpi>0.) fhMassPtGreater3LpiTC->Fill(invMasspKpiLpi);
1802         }
1803         if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1804           if(invMasspiKpKp>0.) fhMassPtGreater3KpTC->Fill(invMasspiKpKp,0.5);
1805           if(invMasspKpiKp>0.) fhMassPtGreater3KpTC->Fill(invMasspKpiKp,0.5);
1806         }else{
1807           if(invMasspiKpKp>0.) fhMassPtGreater3KpTC->Fill(invMasspiKpKp);
1808           if(invMasspKpiKp>0.) fhMassPtGreater3KpTC->Fill(invMasspKpiKp);
1809         }
1810
1811        if(invMasspiKpDk>0. && invMasspKpiDk>0.){
1812         if(invMasspiKpDk>0.) fhMassPtGreater3DkTC->Fill(invMasspiKpDk,0.5);
1813         if(invMasspKpiDk>0.) fhMassPtGreater3DkTC->Fill(invMasspKpiDk,0.5);
1814        }else{
1815         if(invMasspiKpDk>0.) fhMassPtGreater3DkTC->Fill(invMasspiKpDk);
1816         if(invMasspKpiDk>0.) fhMassPtGreater3DkTC->Fill(invMasspKpiDk);
1817        }
1818       if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0.){
1819        if(invMasspiKp3Pr>0.) fhMassPtGreater33PrTC->Fill(invMasspiKp3Pr,0.5);
1820        if(invMasspKpi3Pr>0.) fhMassPtGreater33PrTC->Fill(invMasspKpi3Pr,0.5);
1821       }else{
1822        if(invMasspiKp3Pr>0.) fhMassPtGreater33PrTC->Fill(invMasspiKp3Pr);
1823        if(invMasspKpi3Pr>0.) fhMassPtGreater33PrTC->Fill(invMasspKpi3Pr);
1824       }
1825
1826      }
1827     }
1828     if(part->Pt()>2.&& part->Pt()<=6.){
1829       if(invMasspiKp>0. && invMasspKpi>0.){
1830         if(invMasspiKp>0.) fhMassPtGreater2->Fill(invMasspiKp,0.5);
1831         if(invMasspKpi>0.) fhMassPtGreater2->Fill(invMasspKpi,0.5);
1832       }else{
1833         if(invMasspiKp>0.) fhMassPtGreater2->Fill(invMasspiKp);
1834         if(invMasspKpi>0.) fhMassPtGreater2->Fill(invMasspKpi);
1835       }
1836       if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1837         if(invMasspiKpLpi>0.) fhMassPtGreater2Lpi->Fill(invMasspiKpLpi,0.5);
1838         if(invMasspKpiLpi>0.) fhMassPtGreater2Lpi->Fill(invMasspKpiLpi,0.5);
1839       }else{
1840         if(invMasspiKpLpi>0.) fhMassPtGreater2Lpi->Fill(invMasspiKpLpi);
1841         if(invMasspKpiLpi>0.) fhMassPtGreater2Lpi->Fill(invMasspKpiLpi);
1842       }
1843       if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1844         if(invMasspiKpKp>0.) fhMassPtGreater2Kp->Fill(invMasspiKpKp,0.5);
1845         if(invMasspKpiKp>0.) fhMassPtGreater2Kp->Fill(invMasspKpiKp,0.5);
1846       }else{
1847         if(invMasspiKpKp>0.) fhMassPtGreater2Kp->Fill(invMasspiKpKp);
1848         if(invMasspKpiKp>0.) fhMassPtGreater2Kp->Fill(invMasspKpiKp);
1849       }
1850       if(invMasspiKpDk>0. && invMasspKpiDk>0.){
1851        if(invMasspiKpDk>0.) fhMassPtGreater2Dk->Fill(invMasspiKpDk,0.5);
1852        if(invMasspKpiDk>0.) fhMassPtGreater2Dk->Fill(invMasspKpiDk,0.5);
1853       }else{
1854        if(invMasspiKpDk>0.) fhMassPtGreater2Dk->Fill(invMasspiKpDk);
1855        if(invMasspKpiDk>0.) fhMassPtGreater2Dk->Fill(invMasspKpiDk);
1856       }
1857      if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0.){
1858       if(invMasspiKp3Pr>0.) fhMassPtGreater23Pr->Fill(invMasspiKp3Pr,0.5);
1859       if(invMasspKpi3Pr>0.) fhMassPtGreater23Pr->Fill(invMasspKpi3Pr,0.5);
1860      }else{
1861       if(invMasspiKp3Pr>0.) fhMassPtGreater23Pr->Fill(invMasspiKp3Pr);
1862       if(invMasspKpi3Pr>0.) fhMassPtGreater23Pr->Fill(invMasspKpi3Pr);
1863     }
1864
1865       if(passTightCuts>0){
1866         if(invMasspiKp>0. && invMasspKpi>0.){
1867           if(invMasspiKp>0.) fhMassPtGreater2TC->Fill(invMasspiKp,0.5);
1868           if(invMasspKpi>0.) fhMassPtGreater2TC->Fill(invMasspKpi,0.5);
1869         }else{
1870           if(invMasspiKp>0.) fhMassPtGreater2TC->Fill(invMasspiKp);
1871           if(invMasspKpi>0.) fhMassPtGreater2TC->Fill(invMasspKpi);
1872         }
1873         if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1874           if(invMasspiKpLpi>0.) fhMassPtGreater2LpiTC->Fill(invMasspiKpLpi,0.5);
1875           if(invMasspKpiLpi>0.) fhMassPtGreater2LpiTC->Fill(invMasspKpiLpi,0.5);
1876         }else{
1877           if(invMasspiKpLpi>0.) fhMassPtGreater2LpiTC->Fill(invMasspiKpLpi);
1878           if(invMasspKpiLpi>0.) fhMassPtGreater2LpiTC->Fill(invMasspKpiLpi);
1879         }
1880         if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1881           if(invMasspiKpKp>0.) fhMassPtGreater2KpTC->Fill(invMasspiKpKp,0.5);
1882           if(invMasspKpiKp>0.) fhMassPtGreater2KpTC->Fill(invMasspKpiKp,0.5);
1883         }else{
1884           if(invMasspiKpKp>0.) fhMassPtGreater2KpTC->Fill(invMasspiKpKp);
1885           if(invMasspKpiKp>0.) fhMassPtGreater2KpTC->Fill(invMasspKpiKp);
1886         }
1887        if(invMasspiKpDk>0. && invMasspKpiDk>0.){
1888         if(invMasspiKpDk>0.) fhMassPtGreater2DkTC->Fill(invMasspiKpDk,0.5);
1889         if(invMasspKpiDk>0.) fhMassPtGreater2DkTC->Fill(invMasspKpiDk,0.5);
1890        }else{
1891         if(invMasspiKpDk>0.) fhMassPtGreater2DkTC->Fill(invMasspiKpDk);
1892         if(invMasspKpiDk>0.) fhMassPtGreater2DkTC->Fill(invMasspKpiDk);
1893        }
1894        if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0.){
1895         if(invMasspiKp3Pr>0.) fhMassPtGreater23PrTC->Fill(invMasspiKp3Pr,0.5);
1896         if(invMasspKpi3Pr>0.) fhMassPtGreater23PrTC->Fill(invMasspKpi3Pr,0.5);
1897       }else{
1898        if(invMasspiKp3Pr>0.) fhMassPtGreater23PrTC->Fill(invMasspiKp3Pr);
1899        if(invMasspKpi3Pr>0.) fhMassPtGreater23PrTC->Fill(invMasspKpi3Pr);
1900       }
1901
1902     }
1903    }
1904
1905     if(iPtBin>=0){
1906       index=GetHistoIndex(iPtBin);
1907       if(invMasspiKp>0. && invMasspKpi>0.){
1908         if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp,0.5);
1909         if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi,0.5);
1910       }else{
1911         if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
1912         if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
1913       }
1914       if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1915         if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi,0.5);
1916         if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi,0.5);
1917       }else{
1918         if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi);
1919         if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi);
1920       }
1921       if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1922         if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp,0.5);
1923         if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp,0.5);
1924       }else{
1925         if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp);
1926         if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp);
1927       }
1928       if(invMasspiKpDk>0. && invMasspKpiDk>0.){
1929        if(invMasspiKpDk>0.) fMassHistDk[index]->Fill(invMasspiKpDk,0.5);
1930        if(invMasspKpiDk>0.) fMassHistDk[index]->Fill(invMasspKpiDk,0.5);
1931       }else{
1932        if(invMasspiKpDk>0.) fMassHistDk[index]->Fill(invMasspiKpDk);
1933        if(invMasspKpiDk>0.) fMassHistDk[index]->Fill(invMasspKpiDk);
1934       }
1935       if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0.){
1936        if(invMasspiKp3Pr>0.) fMassHist3Pr[index]->Fill(invMasspiKp3Pr,0.5);
1937        if(invMasspKpi3Pr>0.) fMassHist3Pr[index]->Fill(invMasspKpi3Pr,0.5);
1938       }else{
1939        if(invMasspiKp3Pr>0.) fMassHist3Pr[index]->Fill(invMasspiKp3Pr);
1940        if(invMasspKpi3Pr>0.) fMassHist3Pr[index]->Fill(invMasspKpi3Pr);
1941      }
1942
1943       if(passTightCuts>0){
1944         if(invMasspiKp>0. && invMasspKpi>0. && passTightCuts==3){
1945           if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp,0.5);
1946           if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi,0.5);
1947         }else{
1948           if(invMasspiKp>0. && passTightCuts==2) fMassHistTC[index]->Fill(invMasspiKp);
1949           if(invMasspKpi>0. && passTightCuts==1) fMassHistTC[index]->Fill(invMasspKpi);
1950         }
1951         if(invMasspiKpLpi>0. && invMasspKpiLpi>0. && passTightCuts==3){
1952           if(invMasspiKpLpi>0.) fMassHistLpiTC[index]->Fill(invMasspiKpLpi,0.5);
1953           if(invMasspKpiLpi>0.) fMassHistLpiTC[index]->Fill(invMasspKpiLpi,0.5);
1954         }else{
1955           if(invMasspiKpLpi>0. && passTightCuts==2) fMassHistLpiTC[index]->Fill(invMasspiKpLpi);
1956           if(invMasspKpiLpi>0.&& passTightCuts==1) fMassHistLpiTC[index]->Fill(invMasspKpiLpi);
1957         }
1958         if(invMasspiKpKp>0. && invMasspKpiKp>0. && passTightCuts==3){
1959           if(invMasspiKpKp>0.) fMassHistKpTC[index]->Fill(invMasspiKpKp,0.5);
1960           if(invMasspKpiKp>0.) fMassHistKpTC[index]->Fill(invMasspKpiKp,0.5);
1961         }else{
1962           if(invMasspiKpKp>0. && passTightCuts==2) fMassHistKpTC[index]->Fill(invMasspiKpKp);
1963           if(invMasspKpiKp>0.&& passTightCuts==1) fMassHistKpTC[index]->Fill(invMasspKpiKp);
1964         }
1965        if(invMasspiKpDk>0. && invMasspKpiDk>0. && passTightCuts==3){
1966         if(invMasspiKpDk>0.) fMassHistDkTC[index]->Fill(invMasspiKpDk,0.5);
1967         if(invMasspKpiDk>0.) fMassHistDkTC[index]->Fill(invMasspKpiDk,0.5);
1968        }else{
1969         if(invMasspiKpDk>0. && passTightCuts==2) fMassHistDkTC[index]->Fill(invMasspiKpDk);
1970         if(invMasspKpiDk>0.&& passTightCuts==1) fMassHistDkTC[index]->Fill(invMasspKpiDk);
1971        }
1972        if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0. && passTightCuts==3){
1973         if(invMasspiKp3Pr>0.) fMassHist3PrTC[index]->Fill(invMasspiKp3Pr,0.5);
1974         if(invMasspKpi3Pr>0.) fMassHist3PrTC[index]->Fill(invMasspKpi3Pr,0.5);
1975        }else{
1976         if(invMasspiKp3Pr>0. && passTightCuts==2) fMassHist3PrTC[index]->Fill(invMasspiKp3Pr);
1977         if(invMasspKpi3Pr>0.&& passTightCuts==1) fMassHist3PrTC[index]->Fill(invMasspKpi3Pr);
1978       }
1979
1980       }
1981
1982       if(fReadMC){
1983         if(labDp>0) {
1984           index=GetSignalHistoIndex(iPtBin);
1985           if(invMasspiKp>0. && invMasspKpi>0.){
1986             if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp,0.5);
1987             if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi,0.5);
1988           }else{
1989             if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
1990             if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
1991           }
1992           if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1993             if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi,0.5);
1994             if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi,0.5);
1995           }else{
1996             if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi);
1997             if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi);
1998           }
1999           if(invMasspiKpKp>0. && invMasspKpiKp>0.){
2000             if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp,0.5);
2001             if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp,0.5);
2002           }else{
2003             if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp);
2004             if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp);
2005           }
2006           if(invMasspiKpDk>0. && invMasspKpiDk>0.){
2007             if(invMasspiKpDk>0.) fMassHistDk[index]->Fill(invMasspiKpDk,0.5);
2008             if(invMasspKpiDk>0.) fMassHistDk[index]->Fill(invMasspKpiDk,0.5);
2009           }else{
2010             if(invMasspiKpDk>0.) fMassHistDk[index]->Fill(invMasspiKpDk);
2011             if(invMasspKpiDk>0.) fMassHistDk[index]->Fill(invMasspKpiDk);
2012           }
2013           if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0.){
2014             if(invMasspiKp3Pr>0.) fMassHist3Pr[index]->Fill(invMasspiKp3Pr,0.5);
2015             if(invMasspKpi3Pr>0.) fMassHist3Pr[index]->Fill(invMasspKpi3Pr,0.5);
2016           }else{
2017             if(invMasspiKp3Pr>0.) fMassHist3Pr[index]->Fill(invMasspiKp3Pr);
2018             if(invMasspKpi3Pr>0.) fMassHistDk[index]->Fill(invMasspKpi3Pr);
2019           }
2020
2021           if(passTightCuts>0){
2022             if(invMasspiKp>0. && invMasspKpi>0. && passTightCuts==3){
2023               if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp,0.5);
2024               if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi,0.5);
2025             }else{
2026               if(invMasspiKp>0. && passTightCuts==2) fMassHistTC[index]->Fill(invMasspiKp);
2027               if(invMasspKpi>0.&& passTightCuts==1) fMassHistTC[index]->Fill(invMasspKpi);
2028             }
2029             if(invMasspiKpLpi>0. && invMasspKpiLpi>0. && passTightCuts==3){
2030               if(invMasspiKpLpi>0.) fMassHistLpiTC[index]->Fill(invMasspiKpLpi,0.5);
2031               if(invMasspKpiLpi>0.) fMassHistLpiTC[index]->Fill(invMasspKpiLpi,0.5);
2032             }else{
2033               if(invMasspiKpLpi>0. && passTightCuts==2) fMassHistLpiTC[index]->Fill(invMasspiKpLpi);
2034               if(invMasspKpiLpi>0.&& passTightCuts==1) fMassHistLpiTC[index]->Fill(invMasspKpiLpi);
2035             }
2036             if(invMasspiKpKp>0. && invMasspKpiKp>0. && passTightCuts==3){
2037               if(invMasspiKpKp>0.) fMassHistKpTC[index]->Fill(invMasspiKpKp,0.5);
2038               if(invMasspKpiKp>0.) fMassHistKpTC[index]->Fill(invMasspKpiKp,0.5);
2039             }else{
2040               if(invMasspiKpKp>0. && passTightCuts==2) fMassHistKpTC[index]->Fill(invMasspiKpKp);
2041               if(invMasspKpiKp>0.&& passTightCuts==1) fMassHistKpTC[index]->Fill(invMasspKpiKp);
2042             }
2043             if(invMasspiKpDk>0. && invMasspKpiDk>0. && passTightCuts==3){
2044               if(invMasspiKpDk>0.) fMassHistDkTC[index]->Fill(invMasspiKpDk,0.5);
2045               if(invMasspKpiDk>0.) fMassHistDkTC[index]->Fill(invMasspKpiDk,0.5);
2046             }else{
2047               if(invMasspiKpDk>0. && passTightCuts==2) fMassHistDkTC[index]->Fill(invMasspiKpDk);
2048               if(invMasspKpiDk>0.&& passTightCuts==1) fMassHistDkTC[index]->Fill(invMasspKpiDk);
2049             }
2050             if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0. && passTightCuts==3){
2051               if(invMasspiKp3Pr>0.) fMassHist3PrTC[index]->Fill(invMasspiKp3Pr,0.5);
2052               if(invMasspKpi3Pr>0.) fMassHist3PrTC[index]->Fill(invMasspKpi3Pr,0.5);
2053             }else{
2054               if(invMasspiKp3Pr>0. && passTightCuts==2) fMassHist3PrTC[index]->Fill(invMasspiKp3Pr);
2055               if(invMasspKpi3Pr>0.&& passTightCuts==1) fMassHist3PrTC[index]->Fill(invMasspKpi3Pr);
2056             }
2057           }
2058
2059         }else{
2060           index=GetBackgroundHistoIndex(iPtBin);
2061           if(invMasspiKp>0. && invMasspKpi>0.){
2062             fMassHist[index]->Fill(invMasspiKp,0.5);
2063             fMassHist[index]->Fill(invMasspKpi,0.5);
2064           }else{
2065             if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
2066             if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
2067           }
2068           if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
2069             if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi,0.5);
2070             if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi,0.5);
2071           }else{
2072             if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi);
2073             if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi);
2074           }
2075           if(invMasspiKpKp>0. && invMasspKpiKp>0.){
2076             if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp,0.5);
2077             if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp,0.5);
2078           }else{
2079             if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp);
2080             if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp);
2081           }
2082           if(invMasspiKpDk>0. && invMasspKpiDk>0.){
2083             if(invMasspiKpDk>0.) fMassHistDk[index]->Fill(invMasspiKpDk,0.5);
2084             if(invMasspKpiDk>0.) fMassHistDk[index]->Fill(invMasspKpiDk,0.5);
2085           }else{
2086             if(invMasspiKpDk>0.) fMassHistDk[index]->Fill(invMasspiKpDk);
2087             if(invMasspKpiDk>0.) fMassHistDk[index]->Fill(invMasspKpiDk);
2088           }
2089           if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0.){
2090             if(invMasspiKp3Pr>0.) fMassHist3Pr[index]->Fill(invMasspiKp3Pr,0.5);
2091             if(invMasspKpi3Pr>0.) fMassHist3Pr[index]->Fill(invMasspKpi3Pr,0.5);
2092           }else{
2093             if(invMasspiKp3Pr>0.) fMassHist3Pr[index]->Fill(invMasspiKp3Pr);
2094             if(invMasspKpi3Pr>0.) fMassHistDk[index]->Fill(invMasspKpi3Pr);
2095           }
2096           if(invMasspiKp>0. && invMasspKpi>0. && passTightCuts==3){
2097             fMassHistTC[index]->Fill(invMasspiKp,0.5);
2098             fMassHistTC[index]->Fill(invMasspKpi,0.5);
2099           }else{
2100             if(invMasspiKp>0. && passTightCuts==2) fMassHistTC[index]->Fill(invMasspiKp);
2101             if(invMasspKpi>0. && passTightCuts==1) fMassHistTC[index]->Fill(invMasspKpi);
2102           }
2103           if(invMasspiKpLpi>0. && invMasspKpiLpi>0. && passTightCuts==3){
2104             if(invMasspiKpLpi>0.) fMassHistLpiTC[index]->Fill(invMasspiKpLpi,0.5);
2105             if(invMasspKpiLpi>0.) fMassHistLpiTC[index]->Fill(invMasspKpiLpi,0.5);
2106           }else{
2107             if(invMasspiKpLpi>0. && passTightCuts==2) fMassHistLpiTC[index]->Fill(invMasspiKpLpi);
2108             if(invMasspKpiLpi>0.&& passTightCuts==1) fMassHistLpiTC[index]->Fill(invMasspKpiLpi);
2109           }
2110           if(invMasspiKpKp>0. && invMasspKpiKp>0. && passTightCuts==3){
2111             if(invMasspiKpKp>0.) fMassHistKpTC[index]->Fill(invMasspiKpKp,0.5);
2112             if(invMasspKpiKp>0.) fMassHistKpTC[index]->Fill(invMasspKpiKp,0.5);
2113           }else{
2114             if(invMasspiKpKp>0. && passTightCuts==2) fMassHistKpTC[index]->Fill(invMasspiKpKp);
2115             if(invMasspKpiKp>0.&& passTightCuts==1) fMassHistKpTC[index]->Fill(invMasspKpiKp);
2116           }
2117           if(invMasspiKpDk>0. && invMasspKpiDk>0. && passTightCuts==3){
2118             if(invMasspiKpDk>0.) fMassHistDkTC[index]->Fill(invMasspiKpDk,0.5);
2119             if(invMasspKpiDk>0.) fMassHistDkTC[index]->Fill(invMasspKpiDk,0.5);
2120           }else{
2121             if(invMasspiKpDk>0. && passTightCuts==2) fMassHistDkTC[index]->Fill(invMasspiKpDk);
2122             if(invMasspKpiDk>0.&& passTightCuts==1) fMassHistDkTC[index]->Fill(invMasspKpiDk);
2123            }
2124            if(invMasspiKp3Pr>0. && invMasspKpi3Pr>0. && passTightCuts==3){
2125             if(invMasspiKp3Pr>0.) fMassHist3PrTC[index]->Fill(invMasspiKp3Pr,0.5);
2126             if(invMasspKpi3Pr>0.) fMassHist3PrTC[index]->Fill(invMasspKpi3Pr,0.5);
2127            }else{
2128             if(invMasspiKp3Pr>0. && passTightCuts==2) fMassHist3PrTC[index]->Fill(invMasspiKp3Pr);
2129             if(invMasspKpi3Pr>0.&& passTightCuts==1) fMassHist3PrTC[index]->Fill(invMasspKpi3Pr);
2130           }
2131         }
2132
2133       }
2134     }
2135   }
2136   return;
2137 }
2138 //-----------------------
2139 void AliAnalysisTaskSELambdac::FillVarHists(AliAODRecoDecayHF3Prong *part,
2140                                                TClonesArray *arrMC,
2141                                                AliRDHFCutsLctopKpi *cuts,
2142                                                /* TList *listout,*/
2143                                                AliAODEvent* aod) {
2144   //
2145   // function used in UserExec to fill variable histograms analysing MC
2146   //
2147
2148   TString fillthis="";
2149
2150   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
2151   Double_t invmasscut=0.05;
2152
2153   Int_t pdgDgLctopKpi[3]={2212,321,211};
2154   Int_t lab=-9999;
2155   if(fReadMC)
2156     lab=part->MatchToMC(4122,arrMC,3,pdgDgLctopKpi); //return MC particle label if the array corresponds to a Lc, -1 if not (cf. AliAODRecoDecay.cxx)
2157
2158   Int_t isSelectedPID=cuts->IsSelectedPID(part); // 0 rejected, 1 Lc -> p K- pi+ (K at center because different sign is there),
2159                                                  // 2 Lc -> pi+ K- p (K at center because different sign is there),
2160                                                  // 3 both (it should never happen...)
2161
2162   if (isSelectedPID==0)fNentries->Fill(7);
2163   if (isSelectedPID==1)fNentries->Fill(8);
2164   if (isSelectedPID==2)fNentries->Fill(9);
2165   if (isSelectedPID==3)fNentries->Fill(10);
2166
2167   AliAODTrack *prong0=(AliAODTrack*)part->GetDaughter(0);
2168   AliAODTrack *prong1=(AliAODTrack*)part->GetDaughter(1);
2169   AliAODTrack *prong2=(AliAODTrack*)part->GetDaughter(2);
2170   if (!prong0 || !prong1 || !prong2) {
2171     fNentries->Fill(5);
2172     return;
2173   }
2174
2175   Double_t minvLcpKpi = part->InvMassLcpKpi();
2176   Double_t minvLcpiKp = part->InvMassLcpiKp();
2177
2178   //check pdg of the prongs
2179   Int_t labprong[3]={-1,-1,-1};
2180   if(fReadMC){
2181     labprong[0]=prong0->GetLabel();
2182     labprong[1]=prong1->GetLabel();
2183     labprong[2]=prong2->GetLabel();
2184   }
2185
2186   AliAODMCParticle *mcprong=0;
2187   Int_t pdgProngMC[3]={-1,-1,-1};
2188   if(fReadMC) {
2189     for (Int_t iprong=0;iprong<3;iprong++){
2190       if(labprong[iprong]<0) continue;
2191       mcprong = (AliAODMCParticle*)arrMC->At(labprong[iprong]);
2192       pdgProngMC[iprong]=TMath::Abs(mcprong->GetPdgCode());
2193     }
2194   }
2195
2196   Int_t pdgProngPID[3]={-1,-1,-1};
2197   if(isSelectedPID>0){
2198     pdgProngPID[1]=321;
2199     if(isSelectedPID==1) {pdgProngPID[0]=2212;pdgProngPID[2]=211;}
2200     if(isSelectedPID==2) {pdgProngPID[0]=211;pdgProngPID[2]=2212;}
2201   }
2202
2203   //fill hRealTot ---> PID efficiency denominators
2204   cuts->SetUsePID(kFALSE); //PAOLA
2205   Int_t selectionCandidate=cuts->IsSelected(part,AliRDHFCuts::kCandidate,aod);//PAOLA
2206   if(fReadMC && selectionCandidate>0) { // 3prongs candidate x Lc (no PID)
2207     if(lab>0){ // Signal
2208       for (Int_t iprong=0; iprong<3; iprong++) {
2209         switch (pdgProngMC[iprong]) {
2210         case 2212:
2211           fillthis="hpRealTot";
2212           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2213           break;
2214         case 321:
2215           fillthis="hKRealTot";
2216           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2217           break;
2218         case 211:
2219           fillthis="hpiRealTot";
2220           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2221           break;
2222         default:
2223           break;
2224         }
2225       }
2226     } else { // bkg
2227       for (Int_t iprong=0; iprong<3; iprong++) {
2228         switch (pdgProngMC[iprong]) {
2229         case 2212:
2230           fillthis="hbpRealTot";
2231           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2232           break;
2233         case 321:
2234           fillthis="hbKRealTot";
2235           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2236           break;
2237         case 211:
2238           fillthis="hbpiRealTot";
2239           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2240           break;
2241         default:
2242           break;
2243         }
2244       }
2245     }
2246   }
2247
2248
2249   cuts->SetUsePID(kTRUE); //PAOLA
2250   Int_t selection=cuts->IsSelected(part,AliRDHFCuts::kCandidate,aod);
2251
2252   if ( (lab>0 && fReadMC) ||             // signal X MC
2253        (isSelectedPID>0 && !fReadMC) ) { // signal+bkg X real data
2254
2255     fillthis="hMass";
2256     if ( (fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 4122) ||
2257          (!fReadMC && (isSelectedPID>0 && part->Charge()>0)) ) { // Lc
2258       if ( (pdgProngPID[0]==2212) && (pdgProngPID[1]==321) && (pdgProngPID[2]==211) )
2259         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(minvLcpKpi); // signal+bkg X MC and real data
2260       else if ( (pdgProngPID[0]==211) && (pdgProngPID[1]==321) && (pdgProngPID[2]==2212) )
2261         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(minvLcpiKp); // signal+bkg X MC and real data
2262     }
2263     else if ( (fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == -4122) ||
2264               (!fReadMC && (isSelectedPID>0 && part->Charge()<0)) ) { // anti-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
2271     if (selection>0) { // 3prongs candidate x Lc (yes PID)
2272
2273       Double_t ptmax=0.;
2274
2275       for (Int_t iprong=0; iprong<3; iprong++) {
2276         if (part->PtProng(iprong)>ptmax) ptmax=part->PtProng(iprong);
2277
2278         AliAODTrack *prong = (AliAODTrack*)part->GetDaughter(iprong);
2279         AliAODPid *pidObjtrk = (AliAODPid*)prong->GetDetPid();
2280         AliAODPidHF *pidObj = (AliAODPidHF*)cuts->GetPidHF();
2281         Bool_t hasTOF=pidObj->CheckStatus(prong,"TOF");
2282         Bool_t hasTPC=pidObj->CheckStatus(prong,"TPC");
2283         Double_t tofSignal=0.;
2284         Double_t dedxTPC=0.;
2285         Double_t momTOF=0.;
2286         Double_t momTPC=0.;
2287         if(hasTOF) {
2288           momTOF = prong->P();
2289           tofSignal=pidObjtrk->GetTOFsignal();
2290         }
2291         if(hasTPC) {
2292           momTPC = pidObjtrk->GetTPCmomentum();
2293           dedxTPC=pidObjtrk->GetTPCsignal();
2294         }
2295         switch (pdgProngPID[iprong]) {
2296         case 2212:
2297           fillthis="hpTOFSignal";
2298           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(tofSignal);
2299           fillthis="hpTPCSignal";
2300           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(dedxTPC);
2301           fillthis="hpptProng";
2302           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2303           fillthis="hpd0Prong";
2304           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
2305           fillthis="hpSignalVspTPC";
2306           ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,dedxTPC);
2307           fillthis="hpSignalVspTOF";
2308           ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,tofSignal);
2309           AliPID::EParticleType typep;
2310           typep=AliPID::EParticleType(4);
2311           if(hasTPC) {
2312             //Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typep);
2313             Double_t nsigmap = fPIDResponse->NumberOfSigmasTPC(prong,typep);
2314             fillthis="hpSigmaVspTPC";
2315             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,nsigmap);
2316           }
2317           if(hasTOF){
2318             //Double_t nsigma=fUtilPid->NumberOfSigmasTOF(prong,typep);
2319             Double_t nsigma=fPIDResponse->NumberOfSigmasTOF(prong,typep);
2320             fillthis="hpSigmaVspTOF";
2321             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,nsigma);
2322           }
2323
2324           if (fReadMC) { // fill hpIDTot ---> PID contamination denominator
2325                          //      hIDGood, hnoID ---> PID numerators
2326             fillthis="hpIDTot";
2327             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2328             if(pdgProngMC[iprong]==2212) { // id protons
2329               fillthis="hpIDGood";
2330               ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2331             }
2332             else { // misidentified electrons, muons, pions and kaons
2333               fillthis="hnopIDp";
2334               ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2335             }
2336           }
2337           break;
2338
2339         case 321:
2340           fillthis="hKTOFSignal";
2341           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(tofSignal);
2342           fillthis="hKTPCSignal";
2343           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(dedxTPC);
2344           fillthis="hKptProng";
2345           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2346           fillthis="hKd0Prong";
2347           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
2348           fillthis="hKSignalVspTPC";
2349           ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,dedxTPC);
2350           fillthis="hKSignalVspTOF";
2351           ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,tofSignal);
2352           AliPID::EParticleType typek;
2353           typek=AliPID::EParticleType(3);
2354           if(hasTPC) {
2355             //Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typek);
2356             Double_t nsigmap = fPIDResponse->NumberOfSigmasTPC(prong,typek);
2357             fillthis="hKSigmaVspTPC";
2358             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,nsigmap);
2359           }
2360           if(hasTOF){
2361             //Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typek);
2362             Double_t nsigma=fPIDResponse->NumberOfSigmasTOF(prong,typek); // protection against 'old' AODs...
2363             fillthis="hKSigmaVspTOF";
2364             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,nsigma);
2365           }
2366
2367           if (fReadMC) { // fill hKIDTot ---> PID contamination denominator
2368                          //      hIDGood, hnoID ---> PID numerators
2369             fillthis="hKIDTot";
2370             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2371             if(pdgProngMC[iprong]==321) { // id kaons
2372               fillthis="hKIDGood";
2373               ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2374             }
2375             else { // misidentified electrons, muons, pions and protons
2376               fillthis="hnokIDk";
2377               ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2378             }
2379           }
2380           break;
2381
2382         case 211:
2383           fillthis="hpiTOFSignal";
2384           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(tofSignal);
2385           fillthis="hpiTPCSignal";
2386           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(dedxTPC);
2387           fillthis="hpiptProng";
2388           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2389           fillthis="hpid0Prong";
2390           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
2391           fillthis="hpiSignalVspTPC";
2392           ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,dedxTPC);
2393           fillthis="hpiSignalVspTOF";
2394           ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,tofSignal);
2395           AliPID::EParticleType typepi;
2396           typepi=AliPID::EParticleType(2);
2397           if(hasTPC) {
2398             //Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typepi);
2399             Double_t nsigmap = fPIDResponse->NumberOfSigmasTPC(prong,typepi);
2400             fillthis="hpiSigmaVspTPC";
2401             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,nsigmap);
2402           }
2403           if(hasTOF){
2404             //Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typepi);
2405             Double_t nsigma=fPIDResponse->NumberOfSigmasTOF(prong,typepi); // protection against 'old' AODs...
2406             fillthis="hpiSigmaVspTOF";
2407             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,nsigma);
2408           }
2409
2410           if (fReadMC) { // fill hpiIDTot ---> PID contamination denominator
2411                          //      hIDGood, hnoID ---> PID numerators
2412             fillthis="hpiIDTot";
2413             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2414             if(pdgProngMC[iprong]==211) { // id pions
2415               fillthis="hpiIDGood";
2416               ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2417             }
2418             else { // misidentified electrons, muons, kaons and protons
2419               fillthis="hnopiIDpi";
2420               ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2421             }
2422           }
2423           break;
2424
2425         default:
2426           break;
2427         }
2428
2429       } //end loop on prongs
2430
2431       //now histograms where we don't need to check identity
2432       fillthis = "hDist12toPrim";
2433       ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->GetDist12toPrim());
2434       ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->GetDist23toPrim());
2435       fillthis = "hSigmaVert";
2436       ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->GetSigmaVert());
2437       fillthis = "hDCAs";
2438       Double_t dcas[3];
2439       part->GetDCAs(dcas);
2440       for (Int_t idca=0;idca<3;idca++) ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(dcas[idca]);
2441       fillthis = "hCosPointingAngle";
2442       ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->CosPointingAngle());
2443       fillthis = "hDecayLength";
2444       ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->DecayLength());
2445       Double_t sum2=part->Getd0Prong(0)*part->Getd0Prong(0)+
2446         part->Getd0Prong(1)*part->Getd0Prong(1)+
2447         part->Getd0Prong(2)*part->Getd0Prong(2);
2448       fillthis = "hSum2";
2449       ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(sum2);
2450       fillthis = "hptmax";
2451       ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(ptmax);
2452       fillthis="hLcpt";
2453       ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->Pt());
2454
2455     } // end if (selection)
2456
2457   } else if( lab<=0 && fReadMC) { // bkg x MC
2458
2459     fillthis="hbMass";
2460
2461     if (part->Charge()>0) {      //Lc background
2462       if ( (pdgProngPID[0]==2212) && (pdgProngPID[1]==321) && (pdgProngPID[2]==211) )
2463         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(minvLcpKpi); // bkg X MC
2464       else if ( (pdgProngPID[0]==211) && (pdgProngPID[1]==321) && (pdgProngPID[2]==2212) )
2465         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(minvLcpiKp); // bkg X MC
2466     }
2467     else if (part->Charge()<0){ //anti-Lc background
2468       if ( (pdgProngPID[0]==2212) && (pdgProngPID[1]==321) && (pdgProngPID[2]==211) )
2469         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(minvLcpKpi); // bkg X MC
2470       else if ( (pdgProngPID[0]==211) && (pdgProngPID[1]==321) && (pdgProngPID[2]==2212) )
2471         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(minvLcpiKp); // bkg X MC
2472     }
2473
2474
2475     //apply cut on invariant mass on the pair
2476     if (TMath::Abs(minvLcpKpi-mPDG)<invmasscut || TMath::Abs(minvLcpiKp-mPDG)<invmasscut) {
2477
2478       if (selection>0) { // 3prongs candidate x Lc (yes PID)
2479
2480         Double_t ptmax=0.;
2481         for (Int_t iprong=0; iprong<3; iprong++) {
2482           if(part->PtProng(iprong)>ptmax)ptmax=part->PtProng(iprong);
2483
2484           AliAODTrack *prong = (AliAODTrack*)part->GetDaughter(iprong);
2485           AliAODPid *pidObjtrk = (AliAODPid*)prong->GetDetPid();
2486           AliAODPidHF *pidObj = (AliAODPidHF*)cuts->GetPidHF();
2487           Bool_t hasTOF=pidObj->CheckStatus(prong,"TOF");
2488           Bool_t hasTPC=pidObj->CheckStatus(prong,"TPC");
2489           Double_t tofSignal=0.;
2490           Double_t dedxTPC=0.;
2491           Double_t momTOF=0.;
2492           Double_t momTPC=0.;
2493           if(hasTOF) {
2494             momTOF = prong->P();
2495             tofSignal=pidObjtrk->GetTOFsignal();
2496           }
2497           if(hasTPC) {
2498             momTPC = pidObjtrk->GetTPCmomentum();
2499             dedxTPC=pidObjtrk->GetTPCsignal();
2500           }
2501
2502           switch (pdgProngPID[iprong]) {
2503           case 2212:
2504             fillthis="hbpTOFSignal";
2505             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(tofSignal);
2506             fillthis="hbpTPCSignal";
2507             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(dedxTPC);
2508             fillthis="hbpptProng";
2509             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2510             fillthis="hbpd0Prong";
2511             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
2512             fillthis="hbpSignalVspTPC";
2513             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,dedxTPC);
2514             fillthis="hbpSignalVspTOF";
2515             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,tofSignal);
2516             AliPID::EParticleType typep;
2517             typep=AliPID::EParticleType(4);
2518             if(hasTPC) {
2519               //Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typep);
2520               Double_t nsigmap = fPIDResponse->NumberOfSigmasTPC(prong,typep);
2521               fillthis="hbpSigmaVspTPC";
2522               ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,nsigmap);
2523             }
2524             if(hasTOF){
2525               //Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typep);
2526               Double_t nsigma=fPIDResponse->NumberOfSigmasTOF(prong,typep); // protection against 'old' AODs...
2527               fillthis="hbpSigmaVspTOF";
2528               ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,nsigma);
2529             }
2530             if(fReadMC){ // fill hbpIDTot ---> PID contamination denominator
2531                          //      hbIDGood, hbnoID ---> PID numerators
2532               fillthis="hbpIDTot";
2533               ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2534               if(pdgProngMC[iprong]==2212) { // id protons
2535                 fillthis="hbpIDGood";
2536                 ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2537               }
2538               else { // misidentified electrons, muons, pions and kaons
2539                 fillthis="hbnopIDp";
2540                 ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2541               }
2542             }
2543             break;
2544
2545           case 321:
2546             fillthis="hbKTOFSignal";
2547             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(tofSignal);
2548             fillthis="hbKTPCSignal";
2549             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(dedxTPC);
2550             fillthis="hbKptProng";
2551             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2552             fillthis="hbKd0Prong";
2553             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
2554             fillthis="hbKSignalVspTPC";
2555             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,dedxTPC);
2556             fillthis="hbKSignalVspTOF";
2557             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,tofSignal);
2558             AliPID::EParticleType typek;
2559             typek=AliPID::EParticleType(3);
2560             if(hasTPC) {
2561               //Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typek);
2562               Double_t nsigmap = fPIDResponse->NumberOfSigmasTPC(prong,typek);
2563               fillthis="hbKSigmaVspTPC";
2564               ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,nsigmap);
2565             }
2566             if(hasTOF){
2567               //Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typek);
2568               Double_t nsigma=fPIDResponse->NumberOfSigmasTOF(prong,typek); // protection against 'old' AODs...
2569               fillthis="hbKSigmaVspTOF";
2570               ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,nsigma);
2571             }
2572             if (fReadMC) { // fill hbKIDTot ---> PID contamination denominator
2573                            //      hIDGood, hnoID ---> PID numerators
2574               fillthis="hbKIDTot";
2575               ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2576               if(pdgProngMC[iprong]==321) { // id kaons
2577                 fillthis="hbKIDGood";
2578                 ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2579               }
2580               else { // misidentified electrons, muons, pions and protons
2581                 fillthis="hbnokIDk";
2582                 ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2583               }
2584             }
2585             break;
2586
2587           case 211:
2588             fillthis="hbpiTOFSignal";
2589             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(tofSignal);
2590             fillthis="hbpiTPCSignal";
2591             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(dedxTPC);
2592             fillthis="hbpiptProng";
2593             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2594             fillthis="hbpid0Prong";
2595             ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
2596             fillthis="hbpiSignalVspTPC";
2597             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,dedxTPC);
2598             fillthis="hbpiSignalVspTOF";
2599             ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,tofSignal);
2600             AliPID::EParticleType typepi;
2601             typepi=AliPID::EParticleType(2);
2602             if(hasTPC) {
2603               //Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typepi);
2604               Double_t nsigmap = fPIDResponse->NumberOfSigmasTPC(prong,typepi);
2605               fillthis="hbpiSigmaVspTPC";
2606               ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTPC,nsigmap);
2607             }
2608             if(hasTOF){
2609               //Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typepi);
2610               Double_t nsigma=fPIDResponse->NumberOfSigmasTOF(prong,typepi); // protection against 'old' AODs...
2611               fillthis="hbpiSigmaVspTOF";
2612               ((TH2F*)fOutputMC->FindObject(fillthis))->Fill(momTOF,nsigma);
2613             }
2614             if (fReadMC) { // fill hbpiIDTot ---> PID contamination denominator
2615                            //      hIDGood, hnonIDTot ---> PID numerators
2616               fillthis="hbpiIDTot";
2617               ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2618               if(pdgProngMC[iprong]==211) { // id pions
2619                 fillthis="hbpiIDGood";
2620                 ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2621               }
2622               else { // misidentified electrons, muons, kaons and protons
2623                 fillthis="hbnopiIDpi";
2624                 ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->PtProng(iprong));
2625               }
2626             }
2627             break;
2628
2629           default:
2630             break;
2631           }
2632
2633         } //end loop on prongs
2634
2635         //now histograms where we don't need to check identity
2636         fillthis="hbLcpt";
2637         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->Pt());
2638
2639         fillthis = "hbDist12toPrim";
2640         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->GetDist12toPrim());
2641         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->GetDist23toPrim());
2642         fillthis = "hbSigmaVert";
2643         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->GetSigmaVert());
2644         fillthis = "hbDCAs";
2645         Double_t dcas[3];
2646         part->GetDCAs(dcas);
2647         for (Int_t idca=0;idca<3;idca++)
2648           ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(dcas[idca]);
2649         fillthis = "hbCosPointingAngle";
2650         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->CosPointingAngle());
2651         fillthis = "hbDecayLength";
2652         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(part->DecayLength());
2653         Double_t sum2=part->Getd0Prong(0)*part->Getd0Prong(0)+
2654           part->Getd0Prong(1)*part->Getd0Prong(1)+
2655           part->Getd0Prong(2)*part->Getd0Prong(2);
2656         fillthis = "hbSum2";
2657         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(sum2);
2658         fillthis = "hbptmax";
2659         ((TH1F*)fOutputMC->FindObject(fillthis))->Fill(ptmax);
2660
2661       } // end if (selection)
2662
2663     } // end mass cut
2664
2665   } // end background case
2666
2667   return;
2668
2669 }
2670
2671 //-------------------------
2672 void AliAnalysisTaskSELambdac::FillAPrioriConcentrations(AliAODRecoDecayHF3Prong *part,
2673                                                             AliRDHFCutsLctopKpi *cuts,
2674                                                             AliAODEvent* aod,
2675                                                             TClonesArray *arrMC)
2676 {
2677
2678   // FillAPrioriConcentrations
2679   cuts->SetUsePID(kFALSE); //Annalisa
2680   Int_t isSelected3ProngByLc=cuts->IsSelected(part,AliRDHFCuts::kCandidate,aod);
2681
2682   TString fillthis="";
2683
2684   if(isSelected3ProngByLc>0 && fReadMC) {
2685
2686     for (Int_t ii=0; ii<3; ii++) {
2687       AliAODTrack *prongTest=(AliAODTrack*)part->GetDaughter(ii);
2688       if (!prongTest) continue;
2689       Int_t labprongTest = prongTest->GetLabel();
2690       if(labprongTest<0) continue;
2691       AliAODMCParticle *mcprongTest = (AliAODMCParticle*)arrMC->At(labprongTest);
2692
2693       switch (TMath::Abs(mcprongTest->GetPdgCode())) {
2694       case 11:
2695         fillthis="hElIn3Prong";
2696         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2697         break;
2698       case 13:
2699         fillthis="hMuIn3Prong";
2700         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2701         break;
2702       case 211:
2703         fillthis="hPiIn3Prong";
2704         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2705         break;
2706       case 321:
2707         fillthis="hKaIn3Prong";
2708         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2709         break;
2710       case 2212:
2711         fillthis="hPrIn3Prong";
2712         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2713         break;
2714       default:
2715         break;
2716       }
2717     }
2718   }
2719
2720   cuts->SetUsePID(kTRUE); //Annalisa
2721
2722 }
2723
2724 //-----------------------
2725 Bool_t AliAnalysisTaskSELambdac::Is3ProngFromPDG(AliAODRecoDecayHF3Prong *d,
2726                                                     TClonesArray *arrayMC,
2727                                                     Int_t pdgToBeCompared)
2728 {
2729   //
2730   // Returs kTRUE if at least one of tracks in the current 3prong
2731   // cames from particle with pdgCode=pdgToBeCompared
2732   //
2733
2734   Bool_t localFlag = kFALSE;
2735
2736   for (Int_t ii=0;ii<3;ii++) {
2737     AliAODTrack *daugh = (AliAODTrack*)d->GetDaughter(ii);
2738
2739     Int_t lab = daugh->GetLabel();
2740     if(lab<0) continue;
2741     AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab);
2742     //Int_t partPdg = part->GetPdgCode();
2743     //printf(" -------------------------------------- partLab=%d partPdg=%d ---\n", lab, partPdg);
2744
2745     Int_t motherLabel = part->GetMother();
2746     if(motherLabel<0) continue;
2747
2748     AliAODMCParticle *motherPart = 0;
2749     Int_t motherPdg = 0;
2750     while (!localFlag && motherLabel>=0) {
2751       motherPart = (AliAODMCParticle*)arrayMC->At(motherLabel);
2752       motherPdg = motherPart->GetPdgCode();
2753       //printf(" ------------- motherLab=%d motherPdg=%d ---\n", motherLabel, motherPdg);
2754       if (TMath::Abs(motherPdg)==pdgToBeCompared) {
2755         //printf("-------------------------- trovato quark c/cbar\n");
2756         localFlag = kTRUE;
2757       }
2758       else {
2759         motherLabel = motherPart->GetMother();
2760       }
2761     }
2762
2763   }
2764
2765   return localFlag;
2766
2767 }
2768
2769 //-----------------------
2770 Bool_t AliAnalysisTaskSELambdac::IsTrackFromPDG(const AliAODTrack *daugh,
2771                                                    TClonesArray *arrayMC,
2772                                                    Int_t pdgToBeCompared)
2773 {
2774   //
2775   // Returs kTRUE if at the tracks comes
2776   // from particle with pdgCode=pdgToBeCompared
2777   //
2778
2779   Bool_t localFlag = kFALSE;
2780
2781   Int_t lab = daugh->GetLabel();
2782   if(lab<0) return localFlag;
2783   AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab);
2784   //Int_t partPdg = part->GetPdgCode();
2785   //printf(" -------------------------------------- partLab=%d partPdg=%d ---\n", lab, partPdg);
2786
2787   Int_t motherLabel = part->GetMother();
2788   if(motherLabel<0) return localFlag;
2789
2790   AliAODMCParticle *motherPart = 0;
2791   Int_t motherPdg = 0;
2792   while (!localFlag && motherLabel>=0) {
2793     motherPart = (AliAODMCParticle*)arrayMC->At(motherLabel);
2794     motherPdg = motherPart->GetPdgCode();
2795     //printf(" ------------- motherLab=%d motherPdg=%d ---\n", motherLabel, motherPdg);
2796     if (TMath::Abs(motherPdg)==pdgToBeCompared) {
2797       //printf("-------------------------- trovato quark c/cbar\n");
2798       localFlag = kTRUE;
2799     }
2800     else {
2801       motherLabel = motherPart->GetMother();
2802     }
2803   }
2804
2805   return localFlag;
2806
2807 }
2808
2809
2810
2811
2812 //-----------------------
2813 Bool_t AliAnalysisTaskSELambdac::IsThereAGeneratedLc(TClonesArray *arrayMC)
2814 {
2815   //
2816   // Returs kTRUE if there is a lepton related to the LambdaC
2817   //
2818
2819   Bool_t localFlag = kFALSE;
2820
2821   AliAODMCParticle *searchLc;
2822
2823   for (Int_t iii=0; iii<arrayMC->GetEntries(); iii++) {
2824     searchLc = (AliAODMCParticle*)arrayMC->At(iii);
2825     Int_t searchLcpdg =  searchLc->GetPdgCode();
2826     if (TMath::Abs(searchLcpdg) == 4122){
2827       localFlag = kTRUE;
2828       break;
2829     }
2830   }
2831
2832
2833   return localFlag;
2834
2835 }
2836 //_________________________________
2837 Int_t AliAnalysisTaskSELambdac::NumberPrimaries(const AliAODEvent *aods)
2838 {
2839 ///////////////estimate primaries
2840   Int_t counter=0;
2841
2842
2843   TClonesArray *aodtracks=(TClonesArray *)aods->GetTracks();
2844
2845   // for(Int_t ji=0;ji<aods->GetNTracks();ji++)
2846   for(Int_t ji=0;ji<aodtracks->GetEntriesFast();ji++)
2847     {
2848       AliAODTrack*aodTrack=(AliAODTrack*)aodtracks->UncheckedAt(ji);
2849       if(aodTrack->IsPrimaryCandidate()) counter++;
2850     }
2851   return counter;
2852
2853 }
2854 //_________________________________
2855
2856 void AliAnalysisTaskSELambdac::MultiplicityStudies(AliAODRecoDecayHF3Prong *part,
2857                                                       AliRDHFCutsLctopKpi *cuts,
2858                                                       AliAODEvent* aod,
2859                                                       TClonesArray *arrMC,
2860                                                       Bool_t &flag1,Bool_t &flag2,Bool_t &flag3,
2861                                                       Bool_t &flag4, Bool_t &flag5, Bool_t &flag6)
2862 {
2863
2864
2865 // Multiplicity studies
2866
2867   TString fillthis="";
2868
2869   Int_t pdgDgLctopKpi[3]={2212,321,211};
2870   Int_t lab=-9999;
2871   if(fReadMC)
2872     lab=part->MatchToMC(4122,arrMC,3,pdgDgLctopKpi); //return MC particle label if the array corresponds to a Lc, -1 if not (cf. AliAODRecoDecay.cxx)
2873
2874   cuts->SetUsePID(kFALSE); //Annalisa
2875   Int_t isSelected3ProngByLc=cuts->IsSelected(part,AliRDHFCuts::kCandidate,aod);
2876
2877   if(isSelected3ProngByLc>0 && fReadMC) {
2878     flag1 = kTRUE;
2879     if (lab>0) {
2880       flag3 = kTRUE;
2881     }
2882
2883     Bool_t is3ProngFromJPsi = Is3ProngFromPDG(part,arrMC,443);
2884     if (is3ProngFromJPsi) flag6=is3ProngFromJPsi;
2885
2886     Bool_t is3ProngFromC = Is3ProngFromPDG(part,arrMC,4);
2887     if (is3ProngFromC) flag4=is3ProngFromC;
2888
2889     Bool_t is3ProngFromB = Is3ProngFromPDG(part,arrMC,5);
2890     if (is3ProngFromB) flag5=is3ProngFromB;
2891
2892     for (Int_t ii=0; ii<3; ii++) {
2893       AliAODTrack *prongTest=(AliAODTrack*)part->GetDaughter(ii);
2894       if (!prongTest) continue;
2895       Int_t labprongTest = prongTest->GetLabel();
2896       if(labprongTest<0) continue;
2897       AliAODMCParticle *mcprongTest = (AliAODMCParticle*)arrMC->At(labprongTest);
2898
2899       switch (TMath::Abs(mcprongTest->GetPdgCode())) {
2900       case 11:
2901         fillthis="hElIn3Prong";
2902         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2903         if (IsTrackFromPDG(prongTest,arrMC,443)) {
2904           fillthis="hElIn3Prong6";
2905           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2906         }
2907         break;
2908       case 13:
2909         fillthis="hMuIn3Prong";
2910         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2911         if (IsTrackFromPDG(prongTest,arrMC,443)) {
2912           fillthis="hMuIn3Prong6";
2913           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2914         }
2915         break;
2916       case 211:
2917         fillthis="hPiIn3Prong";
2918         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2919         if (IsTrackFromPDG(prongTest,arrMC,443)) {
2920           fillthis="hPiIn3Prong6";
2921           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2922         }
2923         break;
2924       case 321:
2925         fillthis="hKaIn3Prong";
2926         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2927         if (IsTrackFromPDG(prongTest,arrMC,443)) {
2928           fillthis="hKaIn3Prong6";
2929           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2930         }
2931         break;
2932       case 2212:
2933         fillthis="hPrIn3Prong";
2934         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2935         if (IsTrackFromPDG(prongTest,arrMC,443)) {
2936           fillthis="hPrIn3Prong6";
2937           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2938         }
2939         break;
2940       default:
2941         break;
2942       }
2943
2944       if (lab>0) {
2945         switch (TMath::Abs(mcprongTest->GetPdgCode())) {
2946         case 11:
2947           fillthis="hElIn3Prong1";
2948           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2949           break;
2950         case 13:
2951           fillthis="hMuIn3Prong1";
2952           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2953           break;
2954         case 211:
2955           fillthis="hPiIn3Prong1";
2956           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2957           break;
2958         case 321:
2959           fillthis="hKaIn3Prong1";
2960           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2961           break;
2962         case 2212:
2963           fillthis="hPrIn3Prong1";
2964           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2965           break;
2966         default:
2967           break;
2968         }
2969       } else {
2970         switch (TMath::Abs(mcprongTest->GetPdgCode())) {
2971         case 11:
2972           fillthis="hElIn3Prong2";
2973           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2974           if (IsTrackFromPDG(prongTest,arrMC,4)) {
2975             fillthis="hElIn3Prong3";
2976             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2977           }
2978           if (IsTrackFromPDG(prongTest,arrMC,5)) {
2979             fillthis="hElIn3Prong4";
2980             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2981           }
2982           if (!IsTrackFromPDG(prongTest,arrMC,4) && !IsTrackFromPDG(prongTest,arrMC,5)) {
2983             fillthis="hElIn3Prong5";
2984             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2985           }
2986           break;
2987         case 13:
2988           fillthis="hMuIn3Prong2";
2989           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2990           if (IsTrackFromPDG(prongTest,arrMC,4)) {
2991             fillthis="hMuIn3Prong3";
2992             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2993           }
2994           if (IsTrackFromPDG(prongTest,arrMC,5)) {
2995             fillthis="hMuIn3Prong4";
2996             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
2997           }
2998           if (!IsTrackFromPDG(prongTest,arrMC,4) && !IsTrackFromPDG(prongTest,arrMC,5)) {
2999             fillthis="hMuIn3Prong5";
3000             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3001           }
3002           break;
3003         case 211:
3004           fillthis="hPiIn3Prong2";
3005           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3006           if (IsTrackFromPDG(prongTest,arrMC,4)) {
3007             fillthis="hPiIn3Prong3";
3008             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3009           }
3010           if (IsTrackFromPDG(prongTest,arrMC,5)) {
3011             fillthis="hPiIn3Prong4";
3012             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3013           }
3014           if (!IsTrackFromPDG(prongTest,arrMC,4) && !IsTrackFromPDG(prongTest,arrMC,5)) {
3015             fillthis="hPiIn3Prong5";
3016             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3017           }
3018           break;
3019         case 321:
3020           fillthis="hKaIn3Prong2";
3021           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3022           if (IsTrackFromPDG(prongTest,arrMC,4)) {
3023             fillthis="hKaIn3Prong3";
3024             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3025           }
3026           if (IsTrackFromPDG(prongTest,arrMC,5)) {
3027             fillthis="hKaIn3Prong4";
3028             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3029           }
3030           if (!IsTrackFromPDG(prongTest,arrMC,4) && !IsTrackFromPDG(prongTest,arrMC,5)) {
3031             fillthis="hKaIn3Prong5";
3032             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3033           }
3034           break;
3035         case 2212:
3036           fillthis="hPrIn3Prong2";
3037           ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3038           if (IsTrackFromPDG(prongTest,arrMC,4)) {
3039             fillthis="hPrIn3Prong3";
3040             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3041           }
3042           if (IsTrackFromPDG(prongTest,arrMC,5)) {
3043             fillthis="hPrIn3Prong4";
3044             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3045           }
3046           if (!IsTrackFromPDG(prongTest,arrMC,4) && !IsTrackFromPDG(prongTest,arrMC,5)) {
3047             fillthis="hPrIn3Prong5";
3048             ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3049           }
3050           break;
3051         default:
3052           break;
3053         }
3054       }
3055       /*
3056         if (is3ProngFromC) {
3057         switch (TMath::Abs(mcprongTest->GetPdgCode())) {
3058         case 11:
3059         fillthis="hElIn3Prong3";
3060         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3061         break;
3062         case 13:
3063         fillthis="hMuIn3Prong3";
3064         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3065         break;
3066         case 211:
3067         fillthis="hPiIn3Prong3";
3068         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3069         break;
3070         case 321:
3071         fillthis="hKaIn3Prong3";
3072         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3073         break;
3074         case 2212:
3075         fillthis="hPrIn3Prong3";
3076         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3077         break;
3078         default:
3079         break;
3080         }
3081         } else { // !is3ProngFromC
3082         if (is3ProngFromB) {
3083         switch (TMath::Abs(mcprongTest->GetPdgCode())) {
3084         case 11:
3085         fillthis="hElIn3Prong4";
3086         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3087         break;
3088         case 13:
3089         fillthis="hMuIn3Prong4";
3090         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3091         break;
3092         case 211:
3093         fillthis="hPiIn3Prong4";
3094         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3095         break;
3096         case 321:
3097         fillthis="hKaIn3Prong4";
3098         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3099         break;
3100         case 2212:
3101         fillthis="hPrIn3Prong4";
3102         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3103         break;
3104         default:
3105         break;
3106         }
3107         } else {//!is3ProngFromB
3108         switch (TMath::Abs(mcprongTest->GetPdgCode())) {
3109         case 11:
3110         fillthis="hElIn3Prong5";
3111         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3112         break;
3113         case 13:
3114         fillthis="hMuIn3Prong5";
3115         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3116         break;
3117         case 211:
3118         fillthis="hPiIn3Prong5";
3119         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3120         break;
3121         case 321:
3122         fillthis="hKaIn3Prong5";
3123         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3124         break;
3125         case 2212:
3126         fillthis="hPrIn3Prong5";
3127         ((TH1F*)fAPriori->FindObject(fillthis))->Fill(part->PtProng(ii));
3128         break;
3129         default:
3130         break;
3131         }
3132         }
3133         }
3134       */
3135     }
3136   }
3137
3138   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
3139   Double_t invmasscut=0.05;
3140
3141   Double_t minvLcpKpi = part->InvMassLcpKpi();
3142   Double_t minvLcpiKp = part->InvMassLcpiKp();
3143
3144   cuts->SetUsePID(kTRUE); //Annalisa
3145   Int_t isSelected3ProngByLcPID=cuts->IsSelected(part,AliRDHFCuts::kCandidate,aod);
3146   if (isSelected3ProngByLcPID>0) {
3147     if (TMath::Abs(minvLcpKpi-mPDG)<invmasscut || TMath::Abs(minvLcpiKp-mPDG)<invmasscut) {
3148       flag2 = kTRUE;
3149     }
3150   }
3151
3152
3153 }