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