PWGHF/vertexingHF/AliAnalysisTaskSELambdac.cxx

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