PWG3/vertexingHF/AliAnalysisTaskSELambdac.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 /////////////////////////////////////////////////////////////
  19 //
  20 // AliAnalysisTaskSE for the extraction of signal(e.g Lambdac) of heavy flavor
  21 // decay candidates with the MC truth.
  22 // Authors: r.romita@gsi.de
  23 /////////////////////////////////////////////////////////////
  24
  25 #include <TClonesArray.h>
  26 #include <TNtuple.h>
  27 #include <TCanvas.h>
  28 #include <TList.h>
  29 #include <TString.h>
  30 #include <TH1F.h>
  31 #include <TH2F.h>
  32 #include <TDatabasePDG.h>
  33
  34 #include <AliAnalysisDataSlot.h>
  35 #include <AliAnalysisDataContainer.h>
  36 #include "AliAnalysisManager.h"
  37 #include "AliAODHandler.h"
  38 #include "AliAODEvent.h"
  39 #include "AliAODVertex.h"
  40 #include "AliAODTrack.h"
  41 #include "AliAODMCHeader.h"
  42 #include "AliAODMCParticle.h"
  43 #include "AliAODRecoDecayHF3Prong.h"
  44 #include "AliAnalysisVertexingHF.h"
  45 #include "AliAnalysisTaskSE.h"
  46 #include "AliAnalysisTaskSELambdac.h"
  47 #include "AliKFParticle.h"
  48 #include "AliAODPidHF.h"
  49 #include "AliRDHFCutsLctopKpi.h"
  50 #include "AliRDHFCuts.h"
  51 #include "AliKFVertex.h"
  52 #include "AliESDVertex.h"
  53 #include "AliTOFPIDResponse.h"
  54 #include "AliAODpidUtil.h"
  55 #include "AliAODPid.h"
  56 #include "AliInputEventHandler.h"
  57
  58 ClassImp(AliAnalysisTaskSELambdac)
  59
  60
  61 //________________________________________________________________________
  62 AliAnalysisTaskSELambdac::AliAnalysisTaskSELambdac():
  63 AliAnalysisTaskSE(),
  64 fOutput(0),
  65 fHistNEvents(0),
  66 fhChi2(0),
  67 fhMassPtGreater3(0),
  68 fhMassPtGreater3TC(0),
  69 fhMassPtGreater3Kp(0),
  70 fhMassPtGreater3KpTC(0),
  71 fhMassPtGreater3Lpi(0),
  72 fhMassPtGreater3LpiTC(0),
  73 fhMassPtGreater2(0),
  74 fhMassPtGreater2TC(0),
  75 fhMassPtGreater2Kp(0),
  76 fhMassPtGreater2KpTC(0),
  77 fhMassPtGreater2Lpi(0),
  78 fhMassPtGreater2LpiTC(0),
  79 fNtupleLambdac(0),
  80 fUpmasslimit(2.486),
  81 fLowmasslimit(2.086),
  82 fNPtBins(0),
  83 fRDCutsAnalysis(0),
  84 fRDCutsProduction(0),
  85 fListCuts(0),
  86 fFillNtuple(kFALSE),
  87 fReadMC(kFALSE),
  88 fMCPid(kFALSE),
  89 fRealPid(kFALSE),
  90 fResPid(kTRUE),
  91 fUseKF(kFALSE),
  92 fAnalysis(kFALSE),
  93 fVHF(0),
  94 fFillVarHists(kTRUE),
  95 fNentries(0),
  96 fOutputMC(0),
  97 fUtilPid(0)
  98 {
  99    // Default constructor
 100 }
 101
 102 //________________________________________________________________________
 103 AliAnalysisTaskSELambdac::AliAnalysisTaskSELambdac(const char *name,Bool_t fillNtuple,AliRDHFCutsLctopKpi *lccutsana,AliRDHFCutsLctopKpi *lccutsprod):
 104 AliAnalysisTaskSE(name),
 105 fOutput(0),
 106 fHistNEvents(0),
 107 fhChi2(0),
 108 fhMassPtGreater3(0),
 109 fhMassPtGreater3TC(0),
 110 fhMassPtGreater3Kp(0),
 111 fhMassPtGreater3KpTC(0),
 112 fhMassPtGreater3Lpi(0),
 113 fhMassPtGreater3LpiTC(0),
 114 fhMassPtGreater2(0),
 115 fhMassPtGreater2TC(0),
 116 fhMassPtGreater2Kp(0),
 117 fhMassPtGreater2KpTC(0),
 118 fhMassPtGreater2Lpi(0),
 119 fhMassPtGreater2LpiTC(0),
 120 fNtupleLambdac(0),
 121 fUpmasslimit(2.486),
 122 fLowmasslimit(2.086),
 123 fNPtBins(0),
 124 fRDCutsAnalysis(lccutsana),
 125 fRDCutsProduction(lccutsprod),
 126 fListCuts(0),
 127 fFillNtuple(fillNtuple),
 128 fReadMC(kFALSE),
 129 fMCPid(kFALSE),
 130 fRealPid(kTRUE),
 131 fResPid(kFALSE),
 132 fUseKF(kFALSE),
 133 fAnalysis(kFALSE),
 134 fVHF(0),
 135 fFillVarHists(kTRUE),
 136 fNentries(0),
 137 fOutputMC(0),
 138 fUtilPid(0)
 139 {
 140    SetPtBinLimit(fRDCutsAnalysis->GetNPtBins()+1,fRDCutsAnalysis->GetPtBinLimits());
 141   // Default constructor
 142    // Output slot #1 writes into a TList container
 143   DefineOutput(1,TList::Class());  //My private output
 144   DefineOutput(2,TList::Class());
 145   DefineOutput(3,TList::Class());
 146   DefineOutput(4,TH1F::Class());
 147   if(fFillNtuple){
 148     // Output slot #2 writes into a TNtuple container
 149     DefineOutput(5,TNtuple::Class());  //My private output
 150   }
 151 }
 152
 153 //________________________________________________________________________
 154 AliAnalysisTaskSELambdac::~AliAnalysisTaskSELambdac()
 155 {
 156   // Destructor
 157   if (fOutput) {
 158     delete fOutput;
 159     fOutput = 0;
 160   }
 161   if (fOutputMC) {
 162     delete fOutputMC;
 163     fOutputMC = 0;
 164   }
 165
 166   if (fVHF) {
 167     delete fVHF;
 168     fVHF = 0;
 169   }
 170
 171  if(fRDCutsAnalysis){
 172     delete fRDCutsAnalysis;
 173     fRDCutsAnalysis = 0;
 174  }
 175  if(fRDCutsProduction){
 176     delete fRDCutsProduction;
 177     fRDCutsProduction = 0;
 178  }
 179
 180  if (fListCuts) {
 181    delete fListCuts;
 182    fListCuts = 0;
 183  }
 184  if (fNentries){
 185    delete fNentries;
 186    fNentries = 0;
 187    }
 188  if (fUtilPid){
 189    delete fUtilPid;
 190    fUtilPid = 0;
 191    }
 192 }
 193 //_________________________________________________________________
 194 void  AliAnalysisTaskSELambdac::SetMassLimits(Float_t range){
 195   fUpmasslimit = 2.286+range;
 196   fLowmasslimit = 2.286-range;
 197 }
 198 //_________________________________________________________________
 199 void  AliAnalysisTaskSELambdac::SetMassLimits(Float_t lowlimit, Float_t uplimit){
 200   if(uplimit>lowlimit)
 201     {
 202       fUpmasslimit = lowlimit;
 203       fLowmasslimit = uplimit;
 204     }
 205 }
 206
 207
 208 //________________________________________________________________________
 209 void AliAnalysisTaskSELambdac::SetPtBinLimit(Int_t n, Float_t* lim){
 210   // define pt bins for analysis
 211   if(n>kMaxPtBins){
 212     printf("Max. number of Pt bins = %d\n",kMaxPtBins);
 213     fNPtBins=kMaxPtBins;
 214     fArrayBinLimits[0]=0.;
 215     fArrayBinLimits[1]=2.;
 216     fArrayBinLimits[2]=3.;
 217     fArrayBinLimits[3]=4.;
 218     for(Int_t i=4; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
 219   }else{
 220     fNPtBins=n-1;
 221     fArrayBinLimits[0]=lim[0];
 222     for(Int_t i=1; i<fNPtBins+1; i++)
 223       if(lim[i]>fArrayBinLimits[i-1]){
 224         fArrayBinLimits[i]=lim[i];
 225       }
 226       else {
 227         fArrayBinLimits[i]=fArrayBinLimits[i-1];
 228       }
 229     for(Int_t i=fNPtBins; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
 230   }
 231   if(fDebug > 1){
 232     printf("Number of Pt bins = %d\n",fNPtBins);
 233     for(Int_t i=0; i<fNPtBins; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fArrayBinLimits[i],fArrayBinLimits[i+1]);
 234   }
 235 }
 236 //_________________________________________________________________
 237 Double_t  AliAnalysisTaskSELambdac::GetPtBinLimit(Int_t ibin) const{
 238   if(ibin>fNPtBins)return -1;
 239   return fArrayBinLimits[ibin];
 240 }
 241
 242 //_________________________________________________________________
 243 void AliAnalysisTaskSELambdac::Init()
 244 {
 245   // Initialization
 246
 247   if(fDebug > 1) printf("AnalysisTaskSELambdac::Init() \n");
 248
 249   fListCuts=new TList();
 250   fListCuts->SetOwner();
 251
 252   fListCuts->Add(new AliRDHFCutsLctopKpi(*fRDCutsAnalysis));
 253   fListCuts->Add(new AliRDHFCutsLctopKpi(*fRDCutsProduction));
 254   PostData(2,fListCuts);
 255   return;
 256 }
 257
 258 //________________________________________________________________________
 259 void AliAnalysisTaskSELambdac::UserCreateOutputObjects()
 260 {
 261   // Create the output container
 262   //
 263   if(fDebug > 1) printf("AnalysisTaskSELambdac::UserCreateOutputObjects() \n");
 264
 265   // Several histograms are more conveniently managed in a TList
 266   fOutput = new TList();
 267   fOutput->SetOwner();
 268   fOutput->SetName("OutputHistos");
 269
 270   TString hisname;
 271   Int_t index=0;
 272   Int_t indexLS=0;
 273   for(Int_t i=0;i<fNPtBins;i++){
 274
 275     index=GetHistoIndex(i);
 276     indexLS=GetLSHistoIndex(i);
 277
 278     hisname.Form("hMassPt%d",i);
 279     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 280     fMassHist[index]->Sumw2();
 281     hisname.Form("hMassPt%dTC",i);
 282     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 283     fMassHistTC[index]->Sumw2();
 284
 285     hisname.Form("hMassPtLpi%d",i);
 286     fMassHistLpi[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 287     fMassHistLpi[index]->Sumw2();
 288     hisname.Form("hMassPtLpi%dTC",i);
 289     fMassHistLpiTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 290     fMassHistLpiTC[index]->Sumw2();
 291
 292     hisname.Form("hMassPtKp%d",i);
 293     fMassHistKp[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 294     fMassHistKp[index]->Sumw2();
 295     hisname.Form("hMassPtKp%dTC",i);
 296     fMassHistKpTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 297     fMassHistKpTC[index]->Sumw2();
 298 //signal
 299     index=GetSignalHistoIndex(i);
 300     hisname.Form("hSigPt%d",i);
 301     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 302     fMassHist[index]->Sumw2();
 303     hisname.Form("hSigPt%dTC",i);
 304     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 305     fMassHistTC[index]->Sumw2();
 306     hisname.Form("hSigPtLpi%d",i);
 307     fMassHistLpi[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 308     fMassHistLpi[index]->Sumw2();
 309     hisname.Form("hSigPtLpi%dTC",i);
 310     fMassHistLpiTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 311     fMassHistLpiTC[index]->Sumw2();
 312
 313     hisname.Form("hSigPtKp%d",i);
 314     fMassHistKp[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 315     fMassHistKp[index]->Sumw2();
 316     hisname.Form("hSigPtKp%dTC",i);
 317     fMassHistKpTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 318     fMassHistKpTC[index]->Sumw2();
 319
 320     index=GetBackgroundHistoIndex(i);
 321     hisname.Form("hBkgPt%d",i);
 322     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 323     fMassHist[index]->Sumw2();
 324     hisname.Form("hBkgPt%dTC",i);
 325     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 326     fMassHistTC[index]->Sumw2();
 327     hisname.Form("hBkgPtLpi%d",i);
 328     fMassHistLpi[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 329     fMassHistLpi[index]->Sumw2();
 330     hisname.Form("hBkgPtLpi%dTC",i);
 331     fMassHistLpiTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 332     fMassHistLpiTC[index]->Sumw2();
 333
 334     hisname.Form("hBkgPtKp%d",i);
 335     fMassHistKp[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 336     fMassHistKp[index]->Sumw2();
 337     hisname.Form("hBkgPtKp%dTC",i);
 338     fMassHistKpTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 339     fMassHistKpTC[index]->Sumw2();
 340 }
 341
 342   for(Int_t i=0; i<3*fNPtBins; i++){
 343     fOutput->Add(fMassHist[i]);
 344     fOutput->Add(fMassHistTC[i]);
 345     fOutput->Add(fMassHistLpi[i]);
 346     fOutput->Add(fMassHistLpiTC[i]);
 347     fOutput->Add(fMassHistKp[i]);
 348     fOutput->Add(fMassHistKpTC[i]);
 349   }
 350
 351   fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
 352   fHistNEvents->Sumw2();
 353   fHistNEvents->SetMinimum(0);
 354   fOutput->Add(fHistNEvents);
 355
 356   fhChi2 = new TH1F("fhChi2", "Chi2",100,0.,10.);
 357   fhChi2->Sumw2();
 358   fOutput->Add(fhChi2);
 359
 360   fhMassPtGreater3=new TH1F("fhMassPtGreater3","Pt > 3 GeV/c",100,fLowmasslimit,fUpmasslimit);
 361   fhMassPtGreater3->Sumw2();
 362   fOutput->Add(fhMassPtGreater3);
 363   fhMassPtGreater3TC=new TH1F("fhMassPtGreater3TC","Pt > 3 GeV/c",100,fLowmasslimit,fUpmasslimit);
 364   fhMassPtGreater3TC->Sumw2();
 365   fOutput->Add(fhMassPtGreater3TC);
 366   fhMassPtGreater3Kp=new TH1F("fhMassPtGreater3Kp","Pt > 3 GeV/c",100,fLowmasslimit,fUpmasslimit);
 367   fhMassPtGreater3Kp->Sumw2();
 368   fOutput->Add(fhMassPtGreater3Kp);
 369   fhMassPtGreater3KpTC=new TH1F("fhMassPtGreater3KpTC","Pt > 3 GeV/c",100,fLowmasslimit,fUpmasslimit);
 370   fhMassPtGreater3KpTC->Sumw2();
 371   fOutput->Add(fhMassPtGreater3KpTC);
 372   fhMassPtGreater3Lpi=new TH1F("fhMassPtGreater3Lpi","Pt > 3 GeV/c",100,fLowmasslimit,fUpmasslimit);
 373   fhMassPtGreater3Lpi->Sumw2();
 374   fOutput->Add(fhMassPtGreater3Lpi);
 375   fhMassPtGreater3LpiTC=new TH1F("fhMassPtGreater3LpiTC","Pt > 3 GeV/c",100,fLowmasslimit,fUpmasslimit);
 376   fhMassPtGreater3LpiTC->Sumw2();
 377   fOutput->Add(fhMassPtGreater3LpiTC);
 378   fhMassPtGreater2=new TH1F("fhMassPtGreater2","Pt > 2 GeV/c",100,fLowmasslimit,fUpmasslimit);
 379   fhMassPtGreater2->Sumw2();
 380   fOutput->Add(fhMassPtGreater2);
 381   fhMassPtGreater2TC=new TH1F("fhMassPtGreater2TC","Pt > 2 GeV/c",100,fLowmasslimit,fUpmasslimit);
 382   fhMassPtGreater2TC->Sumw2();
 383   fOutput->Add(fhMassPtGreater2TC);
 384   fhMassPtGreater2Kp=new TH1F("fhMassPtGreater2Kp","Pt > 2 GeV/c",100,fLowmasslimit,fUpmasslimit);
 385   fhMassPtGreater2Kp->Sumw2();
 386   fOutput->Add(fhMassPtGreater2Kp);
 387   fhMassPtGreater2KpTC=new TH1F("fhMassPtGreater2KpTC","Pt > 2 GeV/c",100,fLowmasslimit,fUpmasslimit);
 388   fhMassPtGreater2KpTC->Sumw2();
 389   fOutput->Add(fhMassPtGreater2KpTC);
 390   fhMassPtGreater2Lpi=new TH1F("fhMassPtGreater2Lpi","Pt > 2 GeV/c",100,fLowmasslimit,fUpmasslimit);
 391   fhMassPtGreater2Lpi->Sumw2();
 392   fOutput->Add(fhMassPtGreater2Lpi);
 393   fhMassPtGreater2LpiTC=new TH1F("fhMassPtGreater2LpiTC","Pt > 2 GeV/c",100,fLowmasslimit,fUpmasslimit);
 394   fhMassPtGreater2LpiTC->Sumw2();
 395   fOutput->Add(fhMassPtGreater2LpiTC);
 396
 397   fOutputMC = new TList();
 398   fOutputMC->SetOwner();
 399   fOutputMC->SetName("QAMCHistos");
 400
 401 //  const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
 402
 403   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);
 404
 405   //ROS: qui il bin assignment e' modellato su D0 ma sicuramente iv arie
 406   fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 407   fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
 408   fNentries->GetXaxis()->SetBinLabel(3,"nLcSelected");
 409   fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtxS");
 410   fNentries->GetXaxis()->SetBinLabel(5,"ptbin = -1");
 411   fNentries->GetXaxis()->SetBinLabel(6,"no daughter");
 412   fNentries->GetXaxis()->SetBinLabel(7,"nCandSel(Tr)");
 413   fNentries->GetXaxis()->SetBinLabel(8,"PID=0");
 414   fNentries->GetXaxis()->SetBinLabel(9,"PID=1");
 415   fNentries->GetXaxis()->SetBinLabel(10,"PID=2");
 416   fNentries->GetXaxis()->SetBinLabel(11,"PID=3");
 417   fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
 418
 419   hisname.Form("hMass");
 420   TH1F *hMassInv=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 421   fOutputMC->Add(hMassInv);
 422   hisname.Form("hbMass");
 423   TH1F *hBMassInv=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
 424   fOutputMC->Add(hBMassInv);
 425
 426   // proton specific
 427   hisname.Form("hpTOFSignal");
 428   TH1F *hProtonTOFSignal=new TH1F(hisname.Data(),hisname.Data(),100,12000.,50000.0);
 429   fOutputMC->Add(hProtonTOFSignal);
 430   hisname.Form("hbpTOFSignal");
 431   TH1F *hBProtonTOFSignal=new TH1F(hisname.Data(),hisname.Data(),100,12000.,50000.0);
 432   fOutputMC->Add(hBProtonTOFSignal);
 433
 434   hisname.Form("hpTPCSignal");
 435   TH1F *hProtonTPCSignal=new TH1F(hisname.Data(),hisname.Data(),150,0.,150.0);
 436   fOutputMC->Add(hProtonTPCSignal);
 437   hisname.Form("hbpTPCSignal");
 438   TH1F *hBProtonTPCSignal=new TH1F(hisname.Data(),hisname.Data(),150,0.,150.0);
 439   fOutputMC->Add(hBProtonTPCSignal);
 440
 441   hisname.Form("hpptProng");
 442   TH1F *hProtonPtProng=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.0);
 443   fOutputMC->Add(hProtonPtProng);
 444   hisname.Form("hbpptProng");
 445   TH1F *hBProtonPtProng=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.0);
 446   fOutputMC->Add(hBProtonPtProng);
 447
 448   hisname.Form("hpRealTot");
 449   TH1F *hProtonRealTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 450   fOutputMC->Add(hProtonRealTot);
 451   hisname.Form("hbpRealTot");
 452   TH1F *hBProtonRealTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 453   fOutputMC->Add(hBProtonRealTot);
 454   hisname.Form("hpIDTot");
 455   TH1F *hProtonIDTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 456   fOutputMC->Add(hProtonIDTot);
 457   hisname.Form("hpIDGood");
 458   TH1F *hProtonIDGood=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 459   fOutputMC->Add(hProtonIDGood);
 460   hisname.Form("hbpIDGood");
 461   TH1F *hBProtonIDGood=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 462   fOutputMC->Add(hBProtonIDGood);
 463   hisname.Form("hbpIDTot");
 464   TH1F *hBProtonIDTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 465   fOutputMC->Add(hBProtonIDTot);
 466   hisname.Form("hpnonIDTot");
 467   TH1F *hProtonnonIDTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 468   fOutputMC->Add(hProtonnonIDTot);
 469   hisname.Form("hbpnonIDTot");
 470   TH1F *hBProtonnonIDTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 471   fOutputMC->Add(hBProtonnonIDTot);
 472
 473   hisname.Form("hpd0Prong");
 474   TH1F *hProtond0Prong=new TH1F(hisname.Data(),hisname.Data(),100,-0.1,0.1);
 475   fOutputMC->Add(hProtond0Prong);
 476   hisname.Form("hbpd0Prong");
 477   TH1F *hBProtond0Prong=new TH1F(hisname.Data(),hisname.Data(),100,-0.1,0.1);
 478   fOutputMC->Add(hBProtond0Prong);
 479   hisname.Form("hbpSignalVspTOF");
 480   TH2F *hBpSignalVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,12000.,50000.0);
 481   fOutputMC->Add(hBpSignalVspTOF);
 482   hisname.Form("hbpSignalVspTPC");
 483   TH2F *hBpSignalVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,150,0.,150.0);
 484   fOutputMC->Add(hBpSignalVspTPC);
 485   hisname.Form("hpSignalVspTOF");
 486   TH2F *hpSignalVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,12000.,50000.0);
 487   fOutputMC->Add(hpSignalVspTOF);
 488   hisname.Form("hpSignalVspTPC");
 489   TH2F *hpSignalVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,150,0.,150.0);
 490   fOutputMC->Add(hpSignalVspTPC);
 491
 492   hisname.Form("hpSigmaVspTOF");
 493   TH2F *hpSigmaVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 494   fOutputMC->Add(hpSigmaVspTOF);
 495   hisname.Form("hpSigmaVspTPC");
 496   TH2F *hpSigmaVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 497   fOutputMC->Add(hpSigmaVspTPC);
 498   hisname.Form("hbpSigmaVspTOF");
 499   TH2F *hBpSigmaVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 500   fOutputMC->Add(hBpSigmaVspTOF);
 501   hisname.Form("hbpSigmaVspTPC");
 502   TH2F *hBpSigmaVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 503   fOutputMC->Add(hBpSigmaVspTPC);
 504
 505   //kaon specific
 506
 507   hisname.Form("hKTOFSignal");
 508   TH1F *hKaonTOFSignal=new TH1F(hisname.Data(),hisname.Data(),100,12000.,50000.0);
 509   fOutputMC->Add(hKaonTOFSignal);
 510   hisname.Form("hbKTOFSignal");
 511   TH1F *hBKaonTOFSignal=new TH1F(hisname.Data(),hisname.Data(),100,12000.,50000.0);
 512   fOutputMC->Add(hBKaonTOFSignal);
 513
 514   hisname.Form("hKTPCSignal");
 515   TH1F *hKaonTPCSignal=new TH1F(hisname.Data(),hisname.Data(),150,0.,150.0);
 516   fOutputMC->Add(hKaonTPCSignal);
 517   hisname.Form("hbKTPCSignal");
 518   TH1F *hBKaonTPCSignal=new TH1F(hisname.Data(),hisname.Data(),150,0.,150.0);
 519   fOutputMC->Add(hBKaonTPCSignal);
 520
 521   hisname.Form("hKptProng");
 522   TH1F *hKaonPtProng=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.0);
 523   fOutputMC->Add(hKaonPtProng);
 524   hisname.Form("hbKptProng");
 525   TH1F *hBKaonPtProng=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.0);
 526   fOutputMC->Add(hBKaonPtProng);
 527
 528   hisname.Form("hKRealTot");
 529   TH1F *hKaonRealTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 530   fOutputMC->Add(hKaonRealTot);
 531   hisname.Form("hbKRealTot");
 532   TH1F *hBKaonRealTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 533   fOutputMC->Add(hBKaonRealTot);
 534   hisname.Form("hKIDGood");
 535   TH1F *hKaonIDGood=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 536   fOutputMC->Add(hKaonIDGood);
 537   hisname.Form("hKIDTot");
 538   TH1F *hKaonIDTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 539   fOutputMC->Add(hKaonIDTot);
 540   hisname.Form("hbKIDGood");
 541   TH1F *hBKaonIDGood=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 542   fOutputMC->Add(hBKaonIDGood);
 543   hisname.Form("hbKIDTot");
 544   TH1F *hBKaonIDTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 545   fOutputMC->Add(hBKaonIDTot);
 546   hisname.Form("hKnonIDTot");
 547   TH1F *hKaonnonIDTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 548   fOutputMC->Add(hKaonnonIDTot);
 549   hisname.Form("hbKnonIDTot");
 550   TH1F *hBKaonnonIDTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 551   fOutputMC->Add(hBKaonnonIDTot);
 552
 553
 554   hisname.Form("hKd0Prong");
 555   TH1F *hKaond0Prong=new TH1F(hisname.Data(),hisname.Data(),100,-0.1,0.1);
 556   fOutputMC->Add(hKaond0Prong);
 557   hisname.Form("hbKd0Prong");
 558   TH1F *hBKaond0Prong=new TH1F(hisname.Data(),hisname.Data(),100,-0.1,0.1);
 559   fOutputMC->Add(hBKaond0Prong);
 560   hisname.Form("hbKSignalVspTOF");
 561   TH2F *hbKSignalVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,12000.,50000.0);
 562   fOutputMC->Add(hbKSignalVspTOF);
 563   hisname.Form("hbKSignalVspTPC");
 564   TH2F *hbKSignalVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,150,0.,150.0);
 565   fOutputMC->Add(hbKSignalVspTPC);
 566   hisname.Form("hKSignalVspTOF");
 567   TH2F *hKSignalVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,12000.,50000.0);
 568   fOutputMC->Add(hKSignalVspTOF);
 569   hisname.Form("hKSignalVspTPC");
 570   TH2F *hKSignalVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,150,0.,150.0);
 571   fOutputMC->Add(hKSignalVspTPC);
 572
 573   hisname.Form("hKSigmaVspTOF");
 574   TH2F *hKSigmaVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 575   fOutputMC->Add(hKSigmaVspTOF);
 576
 577   hisname.Form("hKSigmaVspTPC");
 578   TH2F *hKSigmaVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 579   fOutputMC->Add(hKSigmaVspTPC);
 580
 581   hisname.Form("hbKSigmaVspTOF");
 582   TH2F *hBKSigmaVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 583   fOutputMC->Add(hBKSigmaVspTOF);
 584   hisname.Form("hbKSigmaVspTPC");
 585   TH2F *hBKSigmaVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 586   fOutputMC->Add(hBKSigmaVspTPC);
 587
 588   // pion specific
 589   hisname.Form("hpiTOFSignal");
 590   TH1F *hPionTOFSignal=new TH1F(hisname.Data(),hisname.Data(),100,12000.,50000.0);
 591   fOutputMC->Add(hPionTOFSignal);
 592   hisname.Form("hbpiTOFSignal");
 593   TH1F *hBPionTOFSignal=new TH1F(hisname.Data(),hisname.Data(),100,12000.,50000.0);
 594   fOutputMC->Add(hBPionTOFSignal);
 595
 596   hisname.Form("hpiTPCSignal");
 597   TH1F *hPionTPCSignal=new TH1F(hisname.Data(),hisname.Data(),100,30.,100.0);
 598   fOutputMC->Add(hPionTPCSignal);
 599   hisname.Form("hbpiTPCSignal");
 600   TH1F *hBPionTPCSignal=new TH1F(hisname.Data(),hisname.Data(),100,30.,100.0);
 601   fOutputMC->Add(hBPionTPCSignal);
 602
 603   hisname.Form("hpiptProng");
 604   TH1F *hPionPtProng=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.0);
 605   fOutputMC->Add(hPionPtProng);
 606   hisname.Form("hbpiptProng");
 607   TH1F *hBPionPtProng=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.0);
 608   fOutputMC->Add(hBPionPtProng);
 609
 610   hisname.Form("hpiRealTot");
 611   TH1F *hPionRealTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 612   fOutputMC->Add(hPionRealTot);
 613   hisname.Form("hbpiRealTot");
 614   TH1F *hBPionRealTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 615   fOutputMC->Add(hBPionRealTot);
 616   hisname.Form("hpiIDGood");
 617   TH1F *hPionIDGood=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 618   fOutputMC->Add(hPionIDGood);
 619   hisname.Form("hpiIDTot");
 620   TH1F *hPionIDTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 621   fOutputMC->Add(hPionIDTot);
 622   hisname.Form("hbpiIDTot");
 623   TH1F *hBPionIDTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 624   fOutputMC->Add(hBPionIDTot);
 625   hisname.Form("hbpiIDGood");
 626   TH1F *hBPionIDGood=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 627   fOutputMC->Add(hBPionIDGood);
 628   hisname.Form("hpinonIDTot");
 629   TH1F *hPionnonIDTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 630   fOutputMC->Add(hPionnonIDTot);
 631   hisname.Form("hbpinonIDTot");
 632   TH1F *hBPionnonIDTot=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.0);
 633   fOutputMC->Add(hBPionnonIDTot);
 634
 635
 636   hisname.Form("hpid0Prong");
 637   TH1F *hPiond0Prong=new TH1F(hisname.Data(),hisname.Data(),100,-0.1,-0.1);
 638   fOutputMC->Add(hPiond0Prong);
 639   hisname.Form("hbpid0Prong");
 640   TH1F *hBPiond0Prong=new TH1F(hisname.Data(),hisname.Data(),100,-0.1,0.1);
 641   fOutputMC->Add(hBPiond0Prong);
 642   hisname.Form("hbpiSignalVspTOF");
 643   TH2F *hbpiSignalVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,12000.,50000.0);
 644   fOutputMC->Add(hbpiSignalVspTOF);
 645   hisname.Form("hbpiSignalVspTPC");
 646   TH2F *hbpiSignalVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,150,0.,150.0);
 647   fOutputMC->Add(hbpiSignalVspTPC);
 648
 649   hisname.Form("hpiSignalVspTOF");
 650   TH2F *hpiSignalVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,12000.,50000.0);
 651   fOutputMC->Add(hpiSignalVspTOF);
 652   hisname.Form("hpiSignalVspTPC");
 653   TH2F *hpiSignalVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,150,0.,150.0);
 654   fOutputMC->Add(hpiSignalVspTPC);
 655
 656   hisname.Form("hbpiSigmaVspTOF");
 657   TH2F *hBpiSigmaVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 658   fOutputMC->Add(hBpiSigmaVspTOF);
 659
 660   hisname.Form("hbpiSigmaVspTPC");
 661   TH2F *hBpiSigmaVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 662   fOutputMC->Add(hBpiSigmaVspTPC);
 663   hisname.Form("hpiSigmaVspTOF");
 664   TH2F *hpiSigmaVspTOF=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 665   fOutputMC->Add(hpiSigmaVspTOF);
 666
 667   hisname.Form("hpiSigmaVspTPC");
 668   TH2F *hpiSigmaVspTPC=new TH2F(hisname.Data(),hisname.Data(),100,0.,5.0,100,-10.0,10.0);
 669   fOutputMC->Add(hpiSigmaVspTPC);
 670
 671   // other generic
 672   hisname.Form("hLcpt");
 673   TH1F *hLcPt=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.0);
 674   fOutputMC->Add(hLcPt);
 675   hisname.Form("hbLcpt");
 676   TH1F *hBLcPt=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.0);
 677   fOutputMC->Add(hBLcPt);
 678   hisname.Form("hDist12toPrim");
 679   TH1F *hDist12Prim=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.0);
 680   fOutputMC->Add(hDist12Prim);
 681   hisname.Form("hbDist12toPrim");
 682   TH1F *hBDist12Prim=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.0);
 683   fOutputMC->Add(hBDist12Prim);
 684
 685   hisname.Form("hSigmaVert");
 686   TH1F *hSigmaVert=new TH1F(hisname.Data(),hisname.Data(),60,0.,0.06);
 687   fOutputMC->Add(hSigmaVert);
 688   hisname.Form("hbSigmaVert");
 689   TH1F *hBSigmaVert=new TH1F(hisname.Data(),hisname.Data(),60,0.,0.06);
 690   fOutputMC->Add(hBSigmaVert);
 691
 692   hisname.Form("hDCAs");
 693   TH1F *hdcas=new TH1F(hisname.Data(),hisname.Data(),200,0.,0.1);
 694   fOutputMC->Add(hdcas);
 695   hisname.Form("hbDCAs");
 696   TH1F *hBdcas=new TH1F(hisname.Data(),hisname.Data(),200,0.,0.1);
 697   fOutputMC->Add(hBdcas);
 698
 699   hisname.Form("hCosPointingAngle");
 700   TH1F *hCosPointingAngle=new TH1F(hisname.Data(),hisname.Data(),40,0.,1.);
 701   fOutputMC->Add(hCosPointingAngle);
 702   hisname.Form("hbCosPointingAngle");
 703   TH1F *hBCosPointingAngle=new TH1F(hisname.Data(),hisname.Data(),40,0.,1.);
 704   fOutputMC->Add(hBCosPointingAngle);
 705
 706   hisname.Form("hDecayLength");
 707   TH1F *hDecayLength=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 708   fOutputMC->Add(hDecayLength);
 709   hisname.Form("hbDecayLength");
 710   TH1F *hBDecayLength=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 711   fOutputMC->Add(hBDecayLength);
 712
 713   hisname.Form("hSum2");
 714   TH1F *hSum2=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 715   fOutputMC->Add(hSum2);
 716   hisname.Form("hbSum2");
 717   TH1F *hBSum2=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 718   fOutputMC->Add(hBSum2);
 719
 720   hisname.Form("hptmax");
 721   TH1F *hPtMax=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 722   fOutputMC->Add(hPtMax);
 723   hisname.Form("hbptmax");
 724   TH1F *hBPtMax=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 725   fOutputMC->Add(hBPtMax);
 726
 727
 728
 729   if(fRDCutsProduction->GetIsUsePID()){
 730   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 731   AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
 732    AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
 733    fRDCutsProduction->GetPidHF()->SetPidResponse(pidResp);
 734    fRDCutsProduction->GetPidpion()->SetPidResponse(pidResp);
 735    fRDCutsProduction->GetPidprot()->SetPidResponse(pidResp);
 736    fUtilPid=new AliAODpidUtil(pidResp);
 737    fUtilPid=new AliAODpidUtil();
 738    }
 739   if(fRDCutsAnalysis->GetIsUsePID()){
 740    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 741    AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
 742    AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
 743   fRDCutsAnalysis->GetPidHF()->SetPidResponse(pidResp);
 744    fRDCutsAnalysis->GetPidpion()->SetPidResponse(pidResp);
 745    fRDCutsAnalysis->GetPidprot()->SetPidResponse(pidResp);
 746   }
 747
 748
 749   PostData(1,fOutput);
 750   if(fFillVarHists) PostData(3,fOutputMC);
 751   PostData(4,fNentries);
 752   if(fFillNtuple){
 753     //OpenFile(3); // 2 is the slot number of the ntuple
 754
 755     fNtupleLambdac = new TNtuple("fNtupleLambdac","D +","pdg:Px:Py:Pz:PtTrue:VxTrue:VyTrue:VzTrue:Ptpi:PtK:Ptpi2:PtRec:PointingAngle:DecLeng:VxRec:VyRec:VzRec:InvMass:sigvert:d0Pi:d0K:d0Pi2:dca:d0square");
 756   PostData(5,fNtupleLambdac);
 757
 758   }
 759
 760   return;
 761 }
 762
 763 //________________________________________________________________________
 764 void AliAnalysisTaskSELambdac::UserExec(Option_t */*option*/)
 765 {
 766   // Execute analysis for current event:
 767   // heavy flavor candidates association to MC truth
 768
 769    AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 770    //tmp
 771  fHistNEvents->Fill(0); // count event
 772   // Post the data already here
 773
 774   TClonesArray *array3Prong = 0;
 775   TClonesArray *arrayLikeSign =0;
 776   if(!aod && AODEvent() && IsStandardAOD()) {
 777     // In case there is an AOD handler writing a standard AOD, use the AOD
 778     // event in memory rather than the input (ESD) event.
 779     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 780     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 781     // have to taken from the AOD event hold by the AliAODExtension
 782     AliAODHandler* aodHandler = (AliAODHandler*)
 783       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 784     if(aodHandler->GetExtensions()) {
 785       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 786       AliAODEvent *aodFromExt = ext->GetAOD();
 787       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 788       arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
 789     }
 790   } else if(aod) {
 791     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 792     arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
 793   }
 794
 795   if(!array3Prong || !aod) {
 796     printf("AliAnalysisTaskSELambdac::UserExec: Charm3Prong branch not found!\n");
 797     return;
 798   }
 799   if(!arrayLikeSign) {
 800     printf("AliAnalysisTaskSELambdac::UserExec: LikeSign3Prong branch not found!\n");
 801   //  return;
 802   }
 803
 804
 805   TClonesArray *arrayMC=0;
 806   AliAODMCHeader *mcHeader=0;
 807
 808   fNentries->Fill(0);
 809   TString trigclass=aod->GetFiredTriggerClasses();
 810     if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fNentries->Fill(14);
 811     if(!fRDCutsProduction->IsEventSelected(aod)) {
 812      if(fRDCutsProduction->GetWhyRejection()==1) // rejected for pileup
 813       fNentries->Fill(13);
 814       return;
 815     }
 816
 817   // AOD primary vertex
 818   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 819   if(!vtx1) return;
 820
 821
 822   // load MC particles
 823   if(fReadMC){
 824
 825     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 826     if(!arrayMC) {
 827       printf("AliAnalysisTaskSELambdac::UserExec: MC particles branch not found!\n");
 828       return;
 829     }
 830
 831   // load MC header
 832     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 833     if(!mcHeader) {
 834     printf("AliAnalysisTaskSELambdac::UserExec: MC header branch not found!\n");
 835     return;
 836     }
 837   }
 838
 839   Int_t n3Prong = array3Prong->GetEntriesFast();
 840
 841
 842   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 843     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 844
 845
 846     Bool_t unsetvtx=kFALSE;
 847     if(!d->GetOwnPrimaryVtx()){
 848       d->SetOwnPrimaryVtx(vtx1);
 849       unsetvtx=kTRUE;
 850     }
 851
 852 //    if(d->GetSelectionMap()) {if(!d->HasSelectionBit(AliRDHFCuts::kLcCuts)) continue;}
 853 //   if(d->GetSelectionMap()) if(!d->HasSelectionBit(AliRDHFCuts::kLcPID)) continue;
 854
 855     Int_t isSelectedTracks = fRDCutsProduction->IsSelected(d,AliRDHFCuts::kTracks,aod);
 856     if(!isSelectedTracks) continue;
 857
 858     if (fRDCutsProduction->IsInFiducialAcceptance(d->Pt(),d->Y(4122))) fNentries->Fill(6);
 859
 860     Int_t ptbin=fRDCutsProduction->PtBin(d->Pt());
 861     if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
 862
 863     FillMassHists(aod,d,arrayMC,fRDCutsProduction);
 864     if(fFillVarHists) FillVarHists(d,arrayMC,fRDCutsProduction,fOutputMC,aod);
 865
 866       /*
 867       //start OS analysis
 868       if(labDp<0)fHistOSbkg->Fill(d->InvMassDplus());
 869       fHistOS->Fill(d->InvMassDplus());
 870       */
 871
 872     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 873   }
 874
 875   PostData(1,fOutput);
 876   if(fFillVarHists) PostData(3,fOutputMC);
 877   PostData(4,fNentries);
 878
 879   return;
 880 }
 881
 882
 883
 884 //________________________________________________________________________
 885 void AliAnalysisTaskSELambdac::Terminate(Option_t */*option*/)
 886 {
 887   // Terminate analysis
 888   //
 889   if(fDebug > 1) printf("AnalysisTaskSELambdac: Terminate() \n");
 890
 891   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 892   if (!fOutput) {
 893     printf("ERROR: fOutput not available\n");
 894     return;
 895   }
 896  //fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
 897
 898  TString hisname;
 899   if(fFillNtuple){
 900     fNtupleLambdac = dynamic_cast<TNtuple*>(GetOutputData(3));
 901   }
 902
 903
 904  return;
 905 }
 906
 907 //________________________________________________________________________
 908 Int_t AliAnalysisTaskSELambdac::MatchToMCLambdac(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC) const{
 909   // check if the candidate is a Lambdac decaying in pKpi or in the resonant channels
 910   Int_t lambdacLab[3]={0,0,0};
 911   Int_t pdgs[3]={0,0,0};
 912   for(Int_t i=0;i<3;i++){
 913    AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
 914    Int_t lab=daugh->GetLabel();
 915    if(lab<0) return 0;
 916    AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab);
 917    if(!part) continue;
 918    pdgs[i]=part->GetPdgCode();
 919    Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
 920    if(partPdgcode==211 || partPdgcode==321 || partPdgcode==2212){
 921         Int_t motherLabel=part->GetMother();
 922              if(motherLabel<0) return 0;
 923    AliAODMCParticle *motherPart = (AliAODMCParticle*)arrayMC->At(motherLabel);
 924    if(!motherPart) continue;
 925         Int_t motherPdg = TMath::Abs(motherPart->GetPdgCode());
 926    if(motherPdg==4122) {
 927         if(GetLambdacDaugh(motherPart,arrayMC)){lambdacLab[i]=motherLabel;continue;}
 928               }
 929    if(motherPdg==313 || motherPdg==2224 || motherPdg==3124){
 930          Int_t GmotherLabel=motherPart->GetMother();
 931                if(GmotherLabel<0) return 0;
 932          AliAODMCParticle *GmotherPart = (AliAODMCParticle*)arrayMC->At(GmotherLabel);
 933          if(!GmotherPart) continue;
 934                Int_t GmotherPdg = TMath::Abs(GmotherPart->GetPdgCode());
 935          if(GmotherPdg==4122) {
 936                 if(GetLambdacDaugh(GmotherPart,arrayMC)) {lambdacLab[i]=GmotherLabel;continue;}
 937               }
 938         }
 939      }
 940   }
 941
 942  if(lambdacLab[0]==lambdacLab[1] && lambdacLab[1]==lambdacLab[2]) {return lambdacLab[0];}
 943  return 0;
 944
 945 }
 946 //------------------------
 947 Bool_t AliAnalysisTaskSELambdac::GetLambdacDaugh(AliAODMCParticle *part,TClonesArray *arrayMC) const{
 948  // check if the particle is a lambdac and if its decay mode is the correct one
 949  Int_t numberOfLambdac=0;
 950  if(TMath::Abs(part->GetPdgCode())!=4122) return kFALSE;
 951  Int_t daugh_tmp[2];
 952  daugh_tmp[0]=part->GetDaughter(0);
 953  daugh_tmp[1]=part->GetDaughter(1);
 954  Int_t nDaugh = (Int_t)part->GetNDaughters();
 955  if(nDaugh<2) return kFALSE;
 956  if(nDaugh>3) return kFALSE;
 957  AliAODMCParticle* pdaugh1 = (AliAODMCParticle*)arrayMC->At(part->GetDaughter(0));
 958    if(!pdaugh1) {return kFALSE;}
 959  Int_t number1 = TMath::Abs(pdaugh1->GetPdgCode());
 960  AliAODMCParticle* pdaugh2 = (AliAODMCParticle*)arrayMC->At(part->GetDaughter(1));
 961    if(!pdaugh2) {return kFALSE;}
 962   Int_t number2 = TMath::Abs(pdaugh2->GetPdgCode());
 963
 964   if(nDaugh==3){
 965    Int_t thirdDaugh=part->GetDaughter(1)-1;
 966             AliAODMCParticle* pdaugh3 = (AliAODMCParticle*)arrayMC->At(thirdDaugh);
 967   Int_t number3 = TMath::Abs(pdaugh3->GetPdgCode());
 968    if((number1==321 && number2==211 && number3==2212) || (number1==211 && number2==321 && number3==2212) || (number1==211 && number2==2212 && number3==321) || (number1==321 && number2==2212 && number3==211) || (number1==2212 && number2==321 && number3==211) || (number1==2212 && number2==211 && number3==321)) {
 969    numberOfLambdac++;
 970    }
 971   }
 972
 973  if(nDaugh==2){
 974
 975   //Lambda resonant
 976
 977   //Lambda -> p K*0
 978   //
 979   Int_t nfiglieK=0;
 980
 981    if((number1==2212 && number2==313)){
 982      nfiglieK=pdaugh2->GetNDaughters();
 983      if(nfiglieK!=2) return kFALSE;
 984      AliAODMCParticle* pdaughK1 = (AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
 985      AliAODMCParticle* pdaughK2 = (AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
 986      if(!pdaughK1) return kFALSE;
 987      if(!pdaughK2) return kFALSE;
 988      if((TMath::Abs(pdaughK1->GetPdgCode())==211 && TMath::Abs(pdaughK2->GetPdgCode())==321) || (TMath::Abs(pdaughK1->GetPdgCode())==321 && TMath::Abs(pdaughK2->GetPdgCode())==211)) numberOfLambdac++;
 989     }
 990
 991    if((number1==313 && number2==2212)){
 992     nfiglieK=pdaugh1->GetNDaughters();
 993     if(nfiglieK!=2) return kFALSE;
 994     AliAODMCParticle* pdaughK1 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
 995     AliAODMCParticle* pdaughK2 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
 996     if(!pdaughK1) return kFALSE;
 997     if(!pdaughK2) return kFALSE;
 998      if((TMath::Abs(pdaughK1->GetPdgCode())==211 && TMath::Abs(pdaughK2->GetPdgCode())==321) || (TMath::Abs(pdaughK1->GetPdgCode())==321 && TMath::Abs(pdaughK2->GetPdgCode())==211)) numberOfLambdac++;
 999    }
1000
1001    //Lambda -> Delta++ k
1002    Int_t nfiglieDelta=0;
1003    if(number1==321 && number2==2224){
1004     nfiglieDelta=pdaugh2->GetNDaughters();
1005     if(nfiglieDelta!=2) return kFALSE;
1006     AliAODMCParticle *pdaughD1=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
1007     AliAODMCParticle *pdaughD2=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
1008     if(!pdaughD1) return kFALSE;
1009     if(!pdaughD2) return kFALSE;
1010     if((TMath::Abs(pdaughD1->GetPdgCode())==211 && TMath::Abs(pdaughD2->GetPdgCode())==2212) || (TMath::Abs(pdaughD1->GetPdgCode())==2212 && TMath::Abs(pdaughD2->GetPdgCode())==211)) numberOfLambdac++;
1011    }
1012    if(number1==2224 && number2==321){
1013     nfiglieDelta=pdaugh1->GetNDaughters();
1014     if(nfiglieDelta!=2) return kFALSE;
1015     AliAODMCParticle* pdaughD1 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
1016     AliAODMCParticle* pdaughD2 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
1017     if(!pdaughD1) return kFALSE;
1018     if(!pdaughD2) return kFALSE;
1019     if((TMath::Abs(pdaughD1->GetPdgCode())==211 && TMath::Abs(pdaughD2->GetPdgCode())==2212) || (TMath::Abs(pdaughD1->GetPdgCode())==2212 && TMath::Abs(pdaughD2->GetPdgCode())==211)) numberOfLambdac++;
1020    }
1021
1022
1023    //Lambdac -> Lambda(1520) pi
1024    Int_t nfiglieLa=0;
1025    if(number1==3124 && number2==211){
1026     nfiglieLa=pdaugh1->GetNDaughters();
1027     if(nfiglieLa!=2) return kFALSE;
1028     AliAODMCParticle *pdaughL1=(AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
1029     AliAODMCParticle *pdaughL2=(AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
1030     if(!pdaughL1) return kFALSE;
1031     if(!pdaughL2) return kFALSE;
1032     if((TMath::Abs(pdaughL1->GetPdgCode())==321 && TMath::Abs(pdaughL2->GetPdgCode())==2212) || (TMath::Abs(pdaughL1->GetPdgCode())==2212 && TMath::Abs(pdaughL2->GetPdgCode())==321)) numberOfLambdac++;
1033    }
1034    if(number1==211 && number2==3124){
1035     nfiglieLa=pdaugh2->GetNDaughters();
1036     if(nfiglieLa!=2) return kFALSE;
1037     AliAODMCParticle *pdaughL1=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
1038     AliAODMCParticle *pdaughL2=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
1039     if(!pdaughL1) return kFALSE;
1040     if(!pdaughL2) return kFALSE;
1041     if((TMath::Abs(pdaughL1->GetPdgCode())==321 && TMath::Abs(pdaughL2->GetPdgCode())==2212) || (TMath::Abs(pdaughL1->GetPdgCode())==2212 && TMath::Abs(pdaughL2->GetPdgCode())==321)) numberOfLambdac++;
1042
1043    }
1044  }
1045
1046  if(numberOfLambdac>0) {return kTRUE;}
1047          return kFALSE;
1048 }
1049 //-----------------------------
1050 Bool_t AliAnalysisTaskSELambdac::IspKpiMC(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC) const{
1051   // Apply MC PID
1052    Int_t lab[3]={0,0,0},pdgs[3]={0,0,0};
1053    for(Int_t i=0;i<3;i++){
1054     AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
1055     lab[i]=daugh->GetLabel();
1056     if(lab[i]<0) return kFALSE;
1057     AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab[i]);
1058     if(!part) return kFALSE;
1059     pdgs[i]=TMath::Abs(part->GetPdgCode());
1060    }
1061
1062    if(pdgs[0]==2212 && pdgs[1]==321 && pdgs[2]==211) return kTRUE;
1063
1064    return kFALSE;
1065 }
1066 //-----------------------------
1067 Bool_t AliAnalysisTaskSELambdac::IspiKpMC(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC) const{
1068
1069   // Apply MC PID
1070    Int_t lab[3]={0,0,0},pdgs[3]={0,0,0};
1071    for(Int_t i=0;i<3;i++){
1072     AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
1073     lab[i]=daugh->GetLabel();
1074     if(lab[i]<0) return kFALSE;
1075     AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab[i]);
1076     if(!part) return kFALSE;
1077     pdgs[i]=TMath::Abs(part->GetPdgCode());
1078    }
1079
1080    if(pdgs[2]==2212 && pdgs[1]==321 && pdgs[0]==211) {return kTRUE;}
1081
1082    return kFALSE;
1083 }
1084 //--------------------------------------
1085 Bool_t AliAnalysisTaskSELambdac::VertexingKF(AliAODRecoDecayHF3Prong *d,Int_t *pdgs,Double_t field) const{
1086  // apply vertexing KF
1087    Int_t iprongs[3]={0,1,2};
1088    Double_t mass[2]={0.,0.};
1089  //topological constr
1090    AliKFParticle *lambdac=d->ApplyVertexingKF(iprongs,3,pdgs,kTRUE,field,mass);
1091    if(!lambdac) return kFALSE;
1092 //  Double_t probTot=TMath::Prob(Lambdac->GetChi2(),Lambdac->GetNDF());
1093 //  if(probTot<fCutsKF[0]) return kFALSE;
1094   if(lambdac->GetChi2()>fCutsKF[0]) return kFALSE;
1095  //mass constr for K*
1096    Int_t ipRes[2];
1097    Int_t pdgres[2];
1098    mass[0]=0.8961;mass[1]=0.03;
1099    if(TMath::Abs(pdgs[0])==211){
1100     ipRes[0]=0;ipRes[1]=1;
1101     pdgres[0]=pdgs[0];pdgres[1]=321;
1102    }
1103    if(TMath::Abs(pdgs[2])==211){
1104     ipRes[0]=2;ipRes[1]=1;
1105     pdgres[0]=pdgs[2];pdgres[1]=321;
1106    }
1107    AliKFParticle *kappaStar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1108
1109   Double_t probKstar=TMath::Prob(kappaStar->GetChi2(),kappaStar->GetNDF());
1110   if(probKstar>fCutsKF[1]) {
1111     AliAODTrack *esdProng1=(AliAODTrack*)d->GetDaughter(ipRes[0]);
1112     AliAODTrack *esdProng2=(AliAODTrack*)d->GetDaughter(ipRes[1]);
1113     AliKFParticle prong1(*esdProng1,pdgres[0]);
1114     AliKFParticle prong2(*esdProng2,pdgres[1]);
1115     if(kappaStar->GetPt()<fCutsKF[2] && prong1.GetAngle(prong2)>fCutsKF[3]) return kFALSE;
1116   }
1117  //mass constr for Lambda
1118    mass[0]=1.520;mass[1]=0.005;
1119    if(TMath::Abs(pdgs[0])==2212){
1120     ipRes[0]=0;ipRes[1]=1;
1121     pdgres[0]=pdgs[0];pdgres[1]=pdgs[1];
1122    }
1123    if(TMath::Abs(pdgs[2])==2212){
1124     ipRes[0]=2;ipRes[1]=1;
1125     pdgres[0]=pdgs[2];pdgres[1]=pdgs[1];
1126    }
1127    AliKFParticle *lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1128   Double_t probLa=TMath::Prob(lambda1520->GetChi2(),lambda1520->GetNDF());
1129   if(probLa>fCutsKF[4]) {
1130     AliAODTrack *esdProng1=(AliAODTrack*)d->GetDaughter(ipRes[0]);
1131     AliAODTrack *esdProng2=(AliAODTrack*)d->GetDaughter(ipRes[1]);
1132     AliKFParticle prong1(*esdProng1,pdgres[0]);
1133     AliKFParticle prong2(*esdProng2,pdgres[1]);
1134     if(lambda1520->GetPt()<fCutsKF[5] && prong1.GetAngle(prong2)>fCutsKF[6]) return kFALSE;
1135   }
1136  //mass constr for Delta
1137    mass[0]=1.2;mass[1]=0.15;
1138    ipRes[0]=0;ipRes[1]=2;
1139    pdgres[0]=pdgs[0];pdgres[1]=pdgs[2];
1140    AliKFParticle *delta=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1141   Double_t probDelta=TMath::Prob(delta->GetChi2(),delta->GetNDF());
1142   if(probDelta>fCutsKF[7]) {
1143     AliAODTrack *esdProng1=(AliAODTrack*)d->GetDaughter(ipRes[0]);
1144     AliAODTrack *esdProng2=(AliAODTrack*)d->GetDaughter(ipRes[1]);
1145     AliKFParticle prong1(*esdProng1,pdgres[0]);
1146     AliKFParticle prong2(*esdProng2,pdgres[1]);
1147     if(delta->GetPt()<fCutsKF[8] && prong1.GetAngle(prong2)>fCutsKF[9]) return kFALSE;
1148   }
1149  return kTRUE;
1150 }
1151 //-------------------------------------
1152 Int_t AliAnalysisTaskSELambdac::IspiKpResonant(AliAODRecoDecayHF3Prong *d,Double_t field) const{
1153
1154  // apply PID using the resonant channels
1155  //if lambda* -> pk
1156         Double_t mass[2]={1.520,0.005};
1157         Int_t ipRes[2]={1,2};
1158         Int_t pdgres[2]={321,2212};
1159         AliKFParticle *lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1160         Double_t probLa=TMath::Prob(lambda1520->GetChi2(),lambda1520->GetNDF());
1161         if(probLa>0.1) return 1;
1162  //K* -> kpi
1163         mass[0]=0.8961;mass[1]=0.03;
1164         ipRes[0]=0;ipRes[1]=1;
1165         pdgres[0]=211;pdgres[1]=321;
1166         AliKFParticle *kstar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1167         Double_t probKa=TMath::Prob(kstar->GetChi2(),kstar->GetNDF());
1168        if(probKa>0.1) return 2;
1169
1170  return 0;
1171
1172 }
1173 //-------------------------------------
1174 Int_t AliAnalysisTaskSELambdac::IspKpiResonant(AliAODRecoDecayHF3Prong *d,Double_t field) const{
1175
1176  // apply PID using the resonant channels
1177  //if lambda* -> pk
1178         Double_t mass[2]={1.520,0.005};
1179         Int_t ipRes[2]={0,1};
1180         Int_t pdgres[2]={2212,321};
1181         AliKFParticle *lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1182         Double_t probLa=TMath::Prob(lambda1520->GetChi2(),lambda1520->GetNDF());
1183         if(probLa>0.1) return 1;
1184  //K* -> kpi
1185         mass[0]=0.8961;mass[1]=0.03;
1186         ipRes[0]=1;ipRes[1]=2;
1187         pdgres[1]=211;pdgres[0]=321;
1188         AliKFParticle *kstar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1189         Double_t probKa=TMath::Prob(kstar->GetChi2(),kstar->GetNDF());
1190         if(probKa>0.1) return 2;
1191
1192  return 0;
1193
1194 }
1195 //---------------------------
1196 void AliAnalysisTaskSELambdac::FillMassHists(AliAODEvent *aod,AliAODRecoDecayHF3Prong *part, TClonesArray *arrayMC, AliRDHFCutsLctopKpi *cuts){
1197
1198  //if MC PID or no PID, unset PID
1199  if(!fRealPid) cuts->SetUsePID(kFALSE);
1200  Int_t selection=cuts->IsSelected(part,AliRDHFCuts::kCandidate,aod);
1201  if(selection>0){
1202  Int_t iPtBin = -1;
1203  Double_t ptCand = part->Pt();
1204  Int_t index=0;
1205
1206  for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
1207   if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
1208  }
1209
1210
1211  Float_t pdgCode=-2;
1212  Float_t pdgCode1=-1;
1213  Float_t pdgCode2=-1;
1214  Int_t labDp=-1;
1215  Float_t deltaPx=0.;
1216  Float_t deltaPy=0.;
1217  Float_t deltaPz=0.;
1218  Float_t truePt=0.;
1219  Float_t xDecay=0.;
1220  Float_t yDecay=0.;
1221  Float_t zDecay=0.;
1222
1223  if(fReadMC){
1224   labDp = MatchToMCLambdac(part,arrayMC);
1225   if(labDp>0){
1226    AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
1227           AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
1228           AliAODMCParticle *dg1 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(1));
1229           deltaPx=partDp->Px()-part->Px();
1230           deltaPy=partDp->Py()-part->Py();
1231           deltaPz=partDp->Pz()-part->Pz();
1232           truePt=partDp->Pt();
1233           xDecay=dg0->Xv();
1234           yDecay=dg0->Yv();
1235           zDecay=dg0->Zv();
1236           pdgCode=TMath::Abs(partDp->GetPdgCode());
1237           pdgCode1=TMath::Abs(dg0->GetPdgCode());
1238           pdgCode2=TMath::Abs(dg1->GetPdgCode());
1239
1240   }else{
1241    pdgCode=-1;
1242   }
1243  }
1244
1245   Double_t invMasspKpi=-1.;
1246   Double_t invMasspiKp=-1.;
1247   Double_t invMasspKpiLpi=-1.;
1248   Double_t invMasspiKpLpi=-1.;
1249   Double_t invMasspKpiKp=-1.;
1250   Double_t invMasspiKpKp=-1.;
1251   Int_t pdgs[3]={0,0,0};
1252   Double_t field=aod->GetMagneticField();
1253 //apply MC PID
1254   if(fReadMC && fMCPid){
1255
1256   if(IspKpiMC(part,arrayMC)) {
1257    invMasspKpi=part->InvMassLcpKpi();
1258    if(fUseKF){
1259     pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
1260     if(!VertexingKF(part,pdgs,field)) invMasspKpi=-1.;
1261    }
1262   }
1263   if(IspiKpMC(part,arrayMC)) {
1264    invMasspiKp=part->InvMassLcpiKp();
1265    if(fUseKF){
1266     pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
1267     if(!VertexingKF(part,pdgs,field)) invMasspiKp=-1.;
1268    }
1269   }
1270  }
1271
1272  // apply realistic PID
1273   if(fRealPid){
1274    if(selection==1 || selection==3) {
1275     invMasspKpi=part->InvMassLcpKpi();
1276     pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
1277     if(fUseKF){
1278      pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
1279      if(!VertexingKF(part,pdgs,field)) invMasspKpi=-1.;
1280     }
1281    }
1282    if(selection>=2) {
1283     invMasspiKp=part->InvMassLcpiKp();
1284     pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
1285     if(fUseKF){
1286      pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
1287      if(!VertexingKF(part,pdgs,field)) invMasspiKp=-1.;
1288     }
1289    }
1290   }
1291   //apply PID using resonances
1292   if(fResPid && fRealPid){
1293    if((selection==3 || selection==1)) {
1294     if(IspKpiResonant(part,field)==1){
1295      invMasspKpiLpi=part->InvMassLcpKpi();
1296     }
1297     if(IspKpiResonant(part,field)==2){
1298      invMasspKpiKp=part->InvMassLcpKpi();
1299     }
1300    }
1301    if(selection>=2) {
1302     if(IspiKpResonant(part,field)==1){
1303      invMasspiKpLpi=part->InvMassLcpiKp();
1304     }
1305     if(IspiKpResonant(part,field)==2){
1306      invMasspiKpKp=part->InvMassLcpiKp();
1307     }
1308    }
1309   }
1310
1311   // no PID
1312   if(!fResPid && !fRealPid && !fMCPid){
1313    if(selection==2 || selection==3) {
1314     invMasspiKp=part->InvMassLcpiKp();
1315     if(fUseKF){
1316      pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
1317      if(!VertexingKF(part,pdgs,field)) invMasspiKp=-1.;
1318     }
1319    }
1320    if(selection==1 || selection==3){
1321     invMasspKpi=part->InvMassLcpKpi();
1322     if(fUseKF){
1323      pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
1324      if(!VertexingKF(part,pdgs,field)) invMasspKpi=-1.;
1325      }
1326     }
1327    }
1328
1329    if(invMasspiKp<0. && invMasspKpi<0.) return;
1330
1331    Int_t passTightCuts=0;
1332    if(fAnalysis) passTightCuts=fRDCutsAnalysis->IsSelected(part,AliRDHFCuts::kCandidate,aod);
1333
1334
1335
1336    Float_t tmp[24];
1337       if(fFillNtuple){
1338         tmp[0]=pdgCode;
1339         tmp[1]=deltaPx;
1340         tmp[2]=deltaPy;
1341         tmp[3]=deltaPz;
1342         tmp[4]=truePt;
1343         tmp[5]=xDecay;
1344         tmp[6]=yDecay;
1345         tmp[7]=zDecay;
1346         if(pdgCode1==2212) {tmp[8]=part->PtProng(0);}else{tmp[8]=0.;}
1347         if(pdgCode1==211) {tmp[10]=part->PtProng(0);}else{tmp[10]=0.;}
1348         tmp[9]=part->PtProng(1);
1349         if(pdgCode2==211) {tmp[10]=part->PtProng(2);}else{tmp[10]=0.;}
1350         tmp[11]=part->Pt();
1351         tmp[12]=part->CosPointingAngle();
1352         tmp[13]=part->DecayLength();
1353         tmp[14]=part->Xv();
1354         tmp[15]=part->Yv();
1355         tmp[16]=part->Zv();
1356         if(invMasspiKp>0.) tmp[17]=invMasspiKp;
1357         if(invMasspKpi>0.) tmp[17]=invMasspKpi;
1358         tmp[18]=part->GetSigmaVert();
1359         tmp[19]=part->Getd0Prong(0);
1360         tmp[20]=part->Getd0Prong(1);
1361         tmp[21]=part->Getd0Prong(2);
1362         tmp[22]=part->GetDCA();
1363         tmp[23]=part->Prodd0d0();
1364         fNtupleLambdac->Fill(tmp);
1365         PostData(5,fNtupleLambdac);
1366       }
1367
1368     if(part->Pt()>3.){
1369      if(invMasspiKp>0. && invMasspKpi>0.){
1370       if(invMasspiKp>0.) fhMassPtGreater3->Fill(invMasspiKp,0.5);
1371       if(invMasspKpi>0.) fhMassPtGreater3->Fill(invMasspKpi,0.5);
1372      }else{
1373       if(invMasspiKp>0.) fhMassPtGreater3->Fill(invMasspiKp);
1374       if(invMasspKpi>0.) fhMassPtGreater3->Fill(invMasspKpi);
1375      }
1376      if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1377       if(invMasspiKpLpi>0.) fhMassPtGreater3Lpi->Fill(invMasspiKpLpi,0.5);
1378       if(invMasspKpiLpi>0.) fhMassPtGreater3Lpi->Fill(invMasspKpiLpi,0.5);
1379      }else{
1380       if(invMasspiKpLpi>0.) fhMassPtGreater3Lpi->Fill(invMasspiKpLpi);
1381       if(invMasspKpiLpi>0.) fhMassPtGreater3Lpi->Fill(invMasspKpiLpi);
1382      }
1383      if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1384       if(invMasspiKpKp>0.) fhMassPtGreater3Kp->Fill(invMasspiKpKp,0.5);
1385       if(invMasspKpiKp>0.) fhMassPtGreater3Kp->Fill(invMasspKpiKp,0.5);
1386      }else{
1387       if(invMasspiKpKp>0.) fhMassPtGreater3Kp->Fill(invMasspiKpKp);
1388       if(invMasspKpiKp>0.) fhMassPtGreater3Kp->Fill(invMasspKpiKp);
1389      }
1390      if(passTightCuts>0){
1391       if(invMasspiKp>0. && invMasspKpi>0.){
1392        if(invMasspiKp>0.) fhMassPtGreater3TC->Fill(invMasspiKp,0.5);
1393        if(invMasspKpi>0.) fhMassPtGreater3TC->Fill(invMasspKpi,0.5);
1394       }else{
1395        if(invMasspiKp>0.) fhMassPtGreater3TC->Fill(invMasspiKp);
1396        if(invMasspKpi>0.) fhMassPtGreater3TC->Fill(invMasspKpi);
1397       }
1398      if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1399       if(invMasspiKpLpi>0.) fhMassPtGreater3LpiTC->Fill(invMasspiKpLpi,0.5);
1400       if(invMasspKpiLpi>0.) fhMassPtGreater3LpiTC->Fill(invMasspKpiLpi,0.5);
1401      }else{
1402       if(invMasspiKpLpi>0.) fhMassPtGreater3LpiTC->Fill(invMasspiKpLpi);
1403       if(invMasspKpiLpi>0.) fhMassPtGreater3LpiTC->Fill(invMasspKpiLpi);
1404      }
1405      if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1406       if(invMasspiKpKp>0.) fhMassPtGreater3KpTC->Fill(invMasspiKpKp,0.5);
1407       if(invMasspKpiKp>0.) fhMassPtGreater3KpTC->Fill(invMasspKpiKp,0.5);
1408      }else{
1409       if(invMasspiKpKp>0.) fhMassPtGreater3KpTC->Fill(invMasspiKpKp);
1410       if(invMasspKpiKp>0.) fhMassPtGreater3KpTC->Fill(invMasspKpiKp);
1411      }
1412      }
1413     }
1414    if(part->Pt()>2.){
1415      if(invMasspiKp>0. && invMasspKpi>0.){
1416       if(invMasspiKp>0.) fhMassPtGreater2->Fill(invMasspiKp,0.5);
1417       if(invMasspKpi>0.) fhMassPtGreater2->Fill(invMasspKpi,0.5);
1418      }else{
1419       if(invMasspiKp>0.) fhMassPtGreater2->Fill(invMasspiKp);
1420       if(invMasspKpi>0.) fhMassPtGreater2->Fill(invMasspKpi);
1421      }
1422      if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1423       if(invMasspiKpLpi>0.) fhMassPtGreater2Lpi->Fill(invMasspiKpLpi,0.5);
1424       if(invMasspKpiLpi>0.) fhMassPtGreater2Lpi->Fill(invMasspKpiLpi,0.5);
1425      }else{
1426       if(invMasspiKpLpi>0.) fhMassPtGreater2Lpi->Fill(invMasspiKpLpi);
1427       if(invMasspKpiLpi>0.) fhMassPtGreater2Lpi->Fill(invMasspKpiLpi);
1428      }
1429      if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1430       if(invMasspiKpKp>0.) fhMassPtGreater2Kp->Fill(invMasspiKpKp,0.5);
1431       if(invMasspKpiKp>0.) fhMassPtGreater2Kp->Fill(invMasspKpiKp,0.5);
1432      }else{
1433       if(invMasspiKpKp>0.) fhMassPtGreater2Kp->Fill(invMasspiKpKp);
1434       if(invMasspKpiKp>0.) fhMassPtGreater2Kp->Fill(invMasspKpiKp);
1435      }
1436      if(passTightCuts>0){
1437       if(invMasspiKp>0. && invMasspKpi>0.){
1438        if(invMasspiKp>0.) fhMassPtGreater2TC->Fill(invMasspiKp,0.5);
1439        if(invMasspKpi>0.) fhMassPtGreater2TC->Fill(invMasspKpi,0.5);
1440       }else{
1441        if(invMasspiKp>0.) fhMassPtGreater2TC->Fill(invMasspiKp);
1442        if(invMasspKpi>0.) fhMassPtGreater2TC->Fill(invMasspKpi);
1443       }
1444      if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1445       if(invMasspiKpLpi>0.) fhMassPtGreater2LpiTC->Fill(invMasspiKpLpi,0.5);
1446       if(invMasspKpiLpi>0.) fhMassPtGreater2LpiTC->Fill(invMasspKpiLpi,0.5);
1447      }else{
1448       if(invMasspiKpLpi>0.) fhMassPtGreater2LpiTC->Fill(invMasspiKpLpi);
1449       if(invMasspKpiLpi>0.) fhMassPtGreater2LpiTC->Fill(invMasspKpiLpi);
1450      }
1451      if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1452       if(invMasspiKpKp>0.) fhMassPtGreater2KpTC->Fill(invMasspiKpKp,0.5);
1453       if(invMasspKpiKp>0.) fhMassPtGreater2KpTC->Fill(invMasspKpiKp,0.5);
1454      }else{
1455       if(invMasspiKpKp>0.) fhMassPtGreater2KpTC->Fill(invMasspiKpKp);
1456       if(invMasspKpiKp>0.) fhMassPtGreater2KpTC->Fill(invMasspKpiKp);
1457      }
1458      }
1459     }
1460
1461   if(iPtBin>=0){
1462    index=GetHistoIndex(iPtBin);
1463    if(invMasspiKp>0. && invMasspKpi>0.){
1464     if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp,0.5);
1465     if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi,0.5);
1466    }else{
1467     if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
1468     if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
1469    }
1470    if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1471     if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi,0.5);
1472     if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi,0.5);
1473    }else{
1474     if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi);
1475     if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi);
1476    }
1477    if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1478     if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp,0.5);
1479     if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp,0.5);
1480    }else{
1481     if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp);
1482     if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp);
1483    }
1484    if(passTightCuts>0){
1485     if(invMasspiKp>0. && invMasspKpi>0. && passTightCuts==3){
1486      if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp,0.5);
1487      if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi,0.5);
1488     }else{
1489      if(invMasspiKp>0. && passTightCuts==2) fMassHistTC[index]->Fill(invMasspiKp);
1490      if(invMasspKpi>0. && passTightCuts==1) fMassHistTC[index]->Fill(invMasspKpi);
1491     }
1492    if(invMasspiKpLpi>0. && invMasspKpiLpi>0. && passTightCuts==3){
1493     if(invMasspiKpLpi>0.) fMassHistLpiTC[index]->Fill(invMasspiKpLpi,0.5);
1494     if(invMasspKpiLpi>0.) fMassHistLpiTC[index]->Fill(invMasspKpiLpi,0.5);
1495    }else{
1496     if(invMasspiKpLpi>0. && passTightCuts==2) fMassHistLpiTC[index]->Fill(invMasspiKpLpi);
1497     if(invMasspKpiLpi>0.&& passTightCuts==1) fMassHistLpiTC[index]->Fill(invMasspKpiLpi);
1498    }
1499    if(invMasspiKpKp>0. && invMasspKpiKp>0. && passTightCuts==3){
1500     if(invMasspiKpKp>0.) fMassHistKpTC[index]->Fill(invMasspiKpKp,0.5);
1501     if(invMasspKpiKp>0.) fMassHistKpTC[index]->Fill(invMasspKpiKp,0.5);
1502    }else{
1503     if(invMasspiKpKp>0. && passTightCuts==2) fMassHistKpTC[index]->Fill(invMasspiKpKp);
1504     if(invMasspKpiKp>0.&& passTightCuts==1) fMassHistKpTC[index]->Fill(invMasspKpiKp);
1505    }
1506    }
1507    if(fReadMC){
1508     if(labDp>0) {
1509      index=GetSignalHistoIndex(iPtBin);
1510      if(invMasspiKp>0. && invMasspKpi>0.){
1511       if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp,0.5);
1512       if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi,0.5);
1513      }else{
1514       if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
1515       if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
1516      }
1517     if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1518      if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi,0.5);
1519      if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi,0.5);
1520     }else{
1521      if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi);
1522      if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi);
1523     }
1524     if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1525      if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp,0.5);
1526      if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp,0.5);
1527     }else{
1528      if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp);
1529      if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp);
1530     }
1531      if(passTightCuts>0){
1532       if(invMasspiKp>0. && invMasspKpi>0. && passTightCuts==3){
1533        if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp,0.5);
1534        if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi,0.5);
1535       }else{
1536        if(invMasspiKp>0. && passTightCuts==2) fMassHistTC[index]->Fill(invMasspiKp);
1537        if(invMasspKpi>0.&& passTightCuts==1) fMassHistTC[index]->Fill(invMasspKpi);
1538       }
1539     if(invMasspiKpLpi>0. && invMasspKpiLpi>0. && passTightCuts==3){
1540      if(invMasspiKpLpi>0.) fMassHistLpiTC[index]->Fill(invMasspiKpLpi,0.5);
1541      if(invMasspKpiLpi>0.) fMassHistLpiTC[index]->Fill(invMasspKpiLpi,0.5);
1542     }else{
1543      if(invMasspiKpLpi>0. && passTightCuts==2) fMassHistLpiTC[index]->Fill(invMasspiKpLpi);
1544      if(invMasspKpiLpi>0.&& passTightCuts==1) fMassHistLpiTC[index]->Fill(invMasspKpiLpi);
1545     }
1546     if(invMasspiKpKp>0. && invMasspKpiKp>0. && passTightCuts==3){
1547      if(invMasspiKpKp>0.) fMassHistKpTC[index]->Fill(invMasspiKpKp,0.5);
1548      if(invMasspKpiKp>0.) fMassHistKpTC[index]->Fill(invMasspKpiKp,0.5);
1549     }else{
1550      if(invMasspiKpKp>0. && passTightCuts==2) fMassHistKpTC[index]->Fill(invMasspiKpKp);
1551      if(invMasspKpiKp>0.&& passTightCuts==1) fMassHistKpTC[index]->Fill(invMasspKpiKp);
1552     }
1553    }
1554
1555    }else{
1556     index=GetBackgroundHistoIndex(iPtBin);
1557     if(invMasspiKp>0. && invMasspKpi>0.){
1558      fMassHist[index]->Fill(invMasspiKp,0.5);
1559      fMassHist[index]->Fill(invMasspKpi,0.5);
1560     }else{
1561      if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
1562      if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
1563     }
1564     if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1565      if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi,0.5);
1566      if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi,0.5);
1567     }else{
1568      if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi);
1569      if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi);
1570     }
1571     if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1572      if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp,0.5);
1573      if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp,0.5);
1574     }else{
1575      if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp);
1576      if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp);
1577     }
1578     if(invMasspiKp>0. && invMasspKpi>0. && passTightCuts==3){
1579      fMassHistTC[index]->Fill(invMasspiKp,0.5);
1580      fMassHistTC[index]->Fill(invMasspKpi,0.5);
1581     }else{
1582      if(invMasspiKp>0. && passTightCuts==2) fMassHistTC[index]->Fill(invMasspiKp);
1583      if(invMasspKpi>0. && passTightCuts==1) fMassHistTC[index]->Fill(invMasspKpi);
1584     }
1585     if(invMasspiKpLpi>0. && invMasspKpiLpi>0. && passTightCuts==3){
1586      if(invMasspiKpLpi>0.) fMassHistLpiTC[index]->Fill(invMasspiKpLpi,0.5);
1587      if(invMasspKpiLpi>0.) fMassHistLpiTC[index]->Fill(invMasspKpiLpi,0.5);
1588     }else{
1589      if(invMasspiKpLpi>0. && passTightCuts==2) fMassHistLpiTC[index]->Fill(invMasspiKpLpi);
1590      if(invMasspKpiLpi>0.&& passTightCuts==1) fMassHistLpiTC[index]->Fill(invMasspKpiLpi);
1591     }
1592     if(invMasspiKpKp>0. && invMasspKpiKp>0. && passTightCuts==3){
1593      if(invMasspiKpKp>0.) fMassHistKpTC[index]->Fill(invMasspiKpKp,0.5);
1594      if(invMasspKpiKp>0.) fMassHistKpTC[index]->Fill(invMasspKpiKp,0.5);
1595     }else{
1596      if(invMasspiKpKp>0. && passTightCuts==2) fMassHistKpTC[index]->Fill(invMasspiKpKp);
1597      if(invMasspKpiKp>0.&& passTightCuts==1) fMassHistKpTC[index]->Fill(invMasspKpiKp);
1598     }
1599    }
1600
1601   }
1602  }
1603 }
1604  return;
1605 }
1606 //-----------------------
1607 void AliAnalysisTaskSELambdac::FillVarHists(AliAODRecoDecayHF3Prong *part, TClonesArray *arrMC, AliRDHFCutsLctopKpi *cuts, TList *listout,AliAODEvent* aod){
1608   //
1609   // function used in UserExec to fill variable histograms analysing MC
1610   //
1611
1612
1613   Int_t pdgDgLctopKpi[3]={2212,321,211};
1614   Int_t lab=-9999;
1615   if(fReadMC) lab=part->MatchToMC(4122,arrMC,3,pdgDgLctopKpi); //return MC particle label if the array corresponds to a Lc, -1 if not (cf. AliAODRecoDecay.cxx)
1616   Int_t isSelectedPID=cuts->IsSelectedPID(part); //0 rejected,1 Lc -> p K- pi+,2 Lc -> pi+ K- p, (K at center because different sign is there)
1617                                                  //3 both (it should never happen...)
1618
1619   if (isSelectedPID==0)fNentries->Fill(7);
1620   if (isSelectedPID==1)fNentries->Fill(8);
1621   if (isSelectedPID==2)fNentries->Fill(9);
1622   if (isSelectedPID==3)fNentries->Fill(10);
1623
1624
1625   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
1626   Double_t minvLcpKpi = part->InvMassLcpKpi();
1627   Double_t minvLcpiKp = part->InvMassLcpiKp();
1628
1629   Double_t invmasscut=0.05;
1630
1631   TString fillthis="";
1632
1633
1634   AliAODTrack *prong0=(AliAODTrack*)part->GetDaughter(0);
1635   AliAODTrack *prong1=(AliAODTrack*)part->GetDaughter(1);
1636   AliAODTrack *prong2=(AliAODTrack*)part->GetDaughter(2);
1637   if (!prong0 || !prong1 || !prong2) {
1638     fNentries->Fill(5);
1639     return;
1640   }
1641
1642   //check pdg of the prongs
1643   Int_t labprong[3];
1644   if(fReadMC){
1645     labprong[0]=prong0->GetLabel();
1646     labprong[1]=prong1->GetLabel();
1647     labprong[2]=prong2->GetLabel();
1648   }
1649   AliAODMCParticle *mcprong=0;
1650
1651   Int_t pdgProng[3]={0,0,0};
1652   Int_t pdgProngID[3]={0,0,0};
1653     if(fReadMC) {
1654      for (Int_t iprong=0;iprong<3;iprong++){
1655       if(labprong[iprong]<0) continue;
1656       mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
1657       pdgProng[iprong]=TMath::Abs(mcprong->GetPdgCode());
1658      }
1659      if(isSelectedPID>0 && fReadMC){
1660       pdgProngID[1]=321;
1661       if(isSelectedPID==1) {pdgProngID[0]=2212;pdgProngID[2]=211;}
1662       if(isSelectedPID==2) {pdgProngID[0]=211;pdgProngID[2]=2212;}
1663      }
1664     } else {
1665       if (isSelectedPID>0) {
1666         pdgProng[1]=321;
1667         if(isSelectedPID==1) {pdgProng[0]=2212;pdgProng[2]=211;}
1668         if(isSelectedPID==2) {pdgProng[0]=211;pdgProng[2]=2212;}
1669       }
1670     }
1671
1672  Int_t selection=cuts->IsSelected(part,AliRDHFCuts::kCandidate,aod);
1673
1674   if((lab>=0 && fReadMC) || (!fReadMC && (isSelectedPID>0)) ){ //signal (from MC or PID)
1675
1676     fillthis="hMass";
1677
1678      if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 4122)
1679         || (!fReadMC && (isSelectedPID>0 && part->Charge()>0))){    //Lc
1680          if ( (pdgProng[0]==2212) && (pdgProng[1]==321) && (pdgProng[2]==211) ) ((TH1F*)listout->FindObject(fillthis))->Fill(minvLcpKpi);
1681          else if ( (pdgProng[0]==211) && (pdgProng[1]==321) && (pdgProng[2]==2212) ) ((TH1F*)listout->FindObject(fillthis))->Fill(minvLcpiKp);
1682     }
1683     else if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == -4122)
1684           || (!fReadMC && (isSelectedPID>0 && part->Charge()<0))){ //anti-Lc
1685       if ( (pdgProng[0]==2212) && (pdgProng[1]==321) && (pdgProng[2]==211) ) ((TH1F*)listout->FindObject(fillthis))->Fill(minvLcpKpi);
1686       else if ( (pdgProng[0]==211) && (pdgProng[1]==321) && (pdgProng[2]==2212) ) ((TH1F*)listout->FindObject(fillthis))->Fill(minvLcpiKp);
1687     }
1688
1689
1690     //apply cut on invariant mass on the pair
1691     if(selection>0){
1692
1693       Double_t ptmax=0.;
1694       for (Int_t iprong=0; iprong<3; iprong++) {
1695         AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(iprong);
1696         if (fReadMC) {
1697          labprong[iprong]=prong->GetLabel();
1698          if(pdgProngID[iprong]==2212) {
1699           fillthis="hpIDTot";
1700           ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1701          }
1702          if(pdgProngID[iprong]==321) {
1703           fillthis="hKIDTot";
1704           ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1705          }
1706          if(pdgProngID[iprong]==211) {
1707           fillthis="hpiIDTot";
1708           ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1709         }
1710        }
1711         AliAODPid *pidObjtrk=prong->GetDetPid();
1712         AliAODPidHF *pidObj=new AliAODPidHF();
1713         Bool_t hasTOF=pidObj->CheckStatus(prong,"TOF");
1714         Bool_t hasTPC=pidObj->CheckStatus(prong,"TPC");
1715         delete pidObj;
1716         Double_t tofSignal=0.;
1717         Double_t dedxTPC=0.;
1718         Double_t momTOF=0.;
1719         Double_t momTPC=0.;
1720         if(hasTOF) {
1721          momTOF = prong->P();
1722          tofSignal=pidObjtrk->GetTOFsignal();
1723         }
1724         if(hasTPC) {
1725          momTPC = pidObjtrk->GetTPCmomentum();
1726          dedxTPC=pidObjtrk->GetTPCsignal();
1727         }
1728         switch (pdgProng[iprong]) {
1729         case 2212:
1730             fillthis="hpTOFSignal";
1731             ((TH1F*)listout->FindObject(fillthis))->Fill(tofSignal);
1732             fillthis="hpTPCSignal";
1733             ((TH1F*)listout->FindObject(fillthis))->Fill(dedxTPC);
1734             fillthis="hpptProng";
1735             ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1736             fillthis="hpd0Prong";
1737             ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
1738             fillthis="hpSignalVspTPC";
1739             ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,dedxTPC);
1740             fillthis="hpSignalVspTOF";
1741             ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,tofSignal);
1742             ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
1743             AliPID::EParticleType typep;
1744             typep=AliPID::EParticleType(4);
1745             if(hasTPC) {
1746              Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typep);
1747              fillthis="hpSigmaVspTPC";
1748              ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,nsigmap);
1749             }
1750             if(hasTOF){
1751              Double_t nsigma=fUtilPid->NumberOfSigmasTOF(prong,typep);
1752              fillthis="hpSigmaVspTOF";
1753              ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,nsigma);
1754             }
1755             if(fReadMC){
1756               // real protons
1757              fillthis="hpRealTot";
1758              ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1759               // id protons
1760               if(pdgProngID[iprong]==2212) {
1761                fillthis="hpIDGood";
1762                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1763               }
1764               // misidentified kaons, pions
1765               if(pdgProngID[iprong]==321) {
1766                fillthis="hKnonIDTot";
1767                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1768               }
1769               if(pdgProngID[iprong]==211) {
1770                fillthis="hpinonIDTot";
1771                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1772               }
1773             }
1774
1775            break;
1776         case 321:
1777             fillthis="hKTOFSignal";
1778             ((TH1F*)listout->FindObject(fillthis))->Fill(tofSignal);
1779             fillthis="hKTPCSignal";
1780             ((TH1F*)listout->FindObject(fillthis))->Fill(dedxTPC);
1781             fillthis="hKptProng";
1782             ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1783             fillthis="hKd0Prong";
1784             ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
1785             fillthis="hKSignalVspTPC";
1786             ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,dedxTPC);
1787             fillthis="hKSignalVspTOF";
1788             ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,tofSignal);
1789             AliPID::EParticleType typek;
1790             typek=AliPID::EParticleType(3);
1791             if(hasTPC) {
1792              Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typek);
1793              fillthis="hKSigmaVspTPC";
1794              ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,nsigmap);
1795             }
1796             if(hasTOF){
1797              Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typek);
1798              fillthis="hKSigmaVspTOF";
1799              ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,nsigma);
1800             }
1801
1802             if(fReadMC){
1803               // real kaons
1804              fillthis="hKRealTot";
1805              ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1806               // id kaons
1807               if(pdgProngID[iprong]==321) {
1808                fillthis="hKIDGood";
1809                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1810               }
1811               // misidentified protons, pions
1812               if(pdgProngID[iprong]==2212) {
1813                fillthis="hpnonIDTot";
1814                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1815               }
1816               if(pdgProngID[iprong]==211) {
1817                fillthis="hpinonIDTot";
1818                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1819               }
1820             }
1821             break;
1822         case 211:
1823             fillthis="hpiTOFSignal";
1824             ((TH1F*)listout->FindObject(fillthis))->Fill(tofSignal);
1825             fillthis="hpiTPCSignal";
1826             ((TH1F*)listout->FindObject(fillthis))->Fill(dedxTPC);
1827             fillthis="hpiptProng";
1828             ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1829             fillthis="hpid0Prong";
1830             ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
1831             fillthis="hpiSignalVspTPC";
1832             ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,dedxTPC);
1833             fillthis="hpiSignalVspTOF";
1834             ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,tofSignal);
1835             AliPID::EParticleType typepi;
1836             typepi=AliPID::EParticleType(2);
1837             if(hasTPC) {
1838              Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typepi);
1839             fillthis="hpiSigmaVspTPC";
1840             ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,nsigmap);
1841             }
1842             if(hasTOF){
1843              Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typepi);
1844              fillthis="hpiSigmaVspTOF";
1845              ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,nsigma);
1846             }
1847
1848             if(fReadMC){
1849               // real pions
1850              fillthis="hpiRealTot";
1851              ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1852               // id pions
1853               if(pdgProngID[iprong]==211) {
1854                fillthis="hpiIDGood";
1855                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1856               }
1857               // misidentified protons, kaons
1858               if(pdgProngID[iprong]==2212) {
1859                fillthis="hpnonIDTot";
1860                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1861               }
1862               if(pdgProngID[iprong]==321) {
1863                fillthis="hKnonIDTot";
1864                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1865               }
1866             }
1867             break;
1868
1869         default:
1870           break;
1871         }
1872         if(part->PtProng(iprong)>ptmax)ptmax=part->PtProng(iprong);
1873
1874       } //end loop on prongs
1875         //now histograms where we don't need to check identity
1876       fillthis = "hDist12toPrim";
1877       ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDist12toPrim());
1878       ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDist23toPrim());
1879       fillthis = "hSigmaVert";
1880       ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetSigmaVert());
1881       fillthis = "hDCAs";
1882       Double_t dcas[3];
1883       part->GetDCAs(dcas);
1884       for (Int_t idca=0;idca<3;idca++) ((TH1F*)listout->FindObject(fillthis))->Fill(dcas[idca]);
1885       fillthis = "hCosPointingAngle";
1886       ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle());
1887       fillthis = "hDecayLength";
1888       ((TH1F*)listout->FindObject(fillthis))->Fill(part->DecayLength());
1889       Double_t sum2=part->Getd0Prong(0)*part->Getd0Prong(0)+part->Getd0Prong(1)*part->Getd0Prong(1)+part->Getd0Prong(2)*part->Getd0Prong(2);
1890       fillthis = "hSum2";
1891       ((TH1F*)listout->FindObject(fillthis))->Fill(sum2);
1892       fillthis = "hptmax";
1893       ((TH1F*)listout->FindObject(fillthis))->Fill(ptmax);
1894       fillthis="hLcpt";
1895       ((TH1F*)listout->FindObject(fillthis))->Fill(part->Pt());
1896
1897     } //end mass cut
1898
1899
1900   } else if( (lab<0) && fReadMC) { // background **ONLY MC**
1901
1902     fillthis="hbMass";
1903
1904      if (part->Charge()>0){    //Lc background
1905          if ( (pdgProng[0]==2212) && (pdgProng[1]==321) && (pdgProng[2]==211) ) ((TH1F*)listout->FindObject(fillthis))->Fill(minvLcpKpi);
1906          else if ( (pdgProng[0]==211) && (pdgProng[1]==321) && (pdgProng[2]==2212) ) ((TH1F*)listout->FindObject(fillthis))->Fill(minvLcpiKp);
1907      }
1908      else if (part->Charge()<0){ //anti-Lc background
1909       if ( (pdgProng[0]==2212) && (pdgProng[1]==321) && (pdgProng[2]==211) ) ((TH1F*)listout->FindObject(fillthis))->Fill(minvLcpKpi);
1910       else if ( (pdgProng[0]==211) && (pdgProng[1]==321) && (pdgProng[2]==2212) ) ((TH1F*)listout->FindObject(fillthis))->Fill(minvLcpiKp);
1911     }
1912
1913
1914     //apply cut on invariant mass on the pair
1915     if(TMath::Abs(minvLcpKpi-mPDG)<invmasscut || TMath::Abs(minvLcpiKp-mPDG)<invmasscut){
1916     if(selection>0){
1917
1918
1919       Double_t ptmax=0.;
1920
1921       for (Int_t iprong=0; iprong<3; iprong++) {
1922         AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(iprong);
1923         if (fReadMC) {
1924         labprong[iprong]=prong->GetLabel();
1925          if(pdgProngID[iprong]==2212) {
1926           fillthis="hbpIDTot";
1927           ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1928          }
1929          if(pdgProngID[iprong]==321) {
1930           fillthis="hbKIDTot";
1931           ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1932          }
1933          if(pdgProngID[iprong]==211) {
1934           fillthis="hbpiIDTot";
1935           ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1936          }
1937         }
1938         AliAODPid *pidObjtrk=prong->GetDetPid();
1939         AliAODPidHF *pidObj=new AliAODPidHF();
1940         Bool_t hasTOF=pidObj->CheckStatus(prong,"TOF");
1941         Bool_t hasTPC=pidObj->CheckStatus(prong,"TPC");
1942         delete pidObj;
1943         Double_t tofSignal=0.;
1944         Double_t dedxTPC=0.;
1945         Double_t momTOF=0.;
1946         Double_t momTPC=0.;
1947         if(hasTOF) {
1948          momTOF = prong->P();
1949          tofSignal=pidObjtrk->GetTOFsignal();
1950         }
1951         if(hasTPC) {
1952          momTPC = pidObjtrk->GetTPCmomentum();
1953          dedxTPC=pidObjtrk->GetTPCsignal();
1954         }
1955
1956         switch (pdgProng[iprong]) {
1957         case 2212:
1958             fillthis="hbpTOFSignal";
1959             ((TH1F*)listout->FindObject(fillthis))->Fill(tofSignal);
1960             fillthis="hbpTPCSignal";
1961             ((TH1F*)listout->FindObject(fillthis))->Fill(dedxTPC);
1962             fillthis="hbpptProng";
1963             ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1964             fillthis="hbpd0Prong";
1965             ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
1966             fillthis="hbpSignalVspTPC";
1967             ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,dedxTPC);
1968             fillthis="hbpSignalVspTOF";
1969             ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,tofSignal);
1970             AliPID::EParticleType typep;
1971             typep=AliPID::EParticleType(4);
1972             if(hasTPC) {
1973              Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typep);
1974              fillthis="hbpSigmaVspTPC";
1975              ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,nsigmap);
1976             }
1977             if(hasTOF){
1978              Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typep);
1979              fillthis="hbpSigmaVspTOF";
1980              ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,nsigma);
1981             }
1982             if(fReadMC){
1983               // real protons
1984              fillthis="hbpRealTot";
1985              ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1986               // id protons
1987               if(pdgProngID[iprong]==2212) {
1988                fillthis="hbpIDGood";
1989                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1990               }
1991               // misidentified kaons, pions
1992               if(pdgProngID[iprong]==321) {
1993                fillthis="hbKnonIDTot";
1994                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1995               }
1996               if(pdgProngID[iprong]==211) {
1997                fillthis="hbpinonIDTot";
1998                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1999               }
2000             }
2001             break;
2002         case 321:
2003             fillthis="hbKTOFSignal";
2004             ((TH1F*)listout->FindObject(fillthis))->Fill(tofSignal);
2005             fillthis="hbKTPCSignal";
2006             ((TH1F*)listout->FindObject(fillthis))->Fill(dedxTPC);
2007             fillthis="hbKptProng";
2008             ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2009             fillthis="hbKd0Prong";
2010             ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
2011             fillthis="hbKSignalVspTPC";
2012             ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,dedxTPC);
2013             fillthis="hbKSignalVspTOF";
2014             ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,tofSignal);
2015             AliPID::EParticleType typek;
2016             typek=AliPID::EParticleType(3);
2017             if(hasTPC) {
2018              Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typek);
2019              fillthis="hbKSigmaVspTPC";
2020              ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,nsigmap);
2021             }
2022             if(hasTOF){
2023              Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typek);
2024              fillthis="hbKSigmaVspTOF";
2025              ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,nsigma);
2026             }
2027             if(fReadMC){
2028               // real kaons
2029              fillthis="hbKRealTot";
2030              ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2031               // id kaons
2032               if(pdgProngID[iprong]==321) {
2033                fillthis="hbKIDGood";
2034                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2035               }
2036               // misidentified protons, pions
2037               if(pdgProngID[iprong]==2212) {
2038                fillthis="hbpnonIDTot";
2039                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2040               }
2041               if(pdgProngID[iprong]==211) {
2042                fillthis="hbpinonIDTot";
2043                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2044               }
2045             }
2046             break;
2047         case 211:
2048             fillthis="hbpiTOFSignal";
2049             ((TH1F*)listout->FindObject(fillthis))->Fill(tofSignal);
2050             fillthis="hbpiTPCSignal";
2051             ((TH1F*)listout->FindObject(fillthis))->Fill(dedxTPC);
2052             fillthis="hbpiptProng";
2053             ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2054             fillthis="hbpid0Prong";
2055             ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
2056             fillthis="hbpiSignalVspTPC";
2057             ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,dedxTPC);
2058             fillthis="hbpiSignalVspTOF";
2059             ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,tofSignal);
2060             AliPID::EParticleType typepi;
2061             typepi=AliPID::EParticleType(2);
2062             if(hasTPC) {
2063              Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typepi);
2064              fillthis="hbpiSigmaVspTPC";
2065              ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,nsigmap);
2066             }
2067             if(hasTOF){
2068              Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typepi);
2069              fillthis="hbpiSigmaVspTOF";
2070              ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,nsigma);
2071             }
2072             if(fReadMC){
2073               // real pions
2074              fillthis="hbpiRealTot";
2075              ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2076               // id pions
2077               if(pdgProngID[iprong]==211) {
2078                fillthis="hbpiIDGood";
2079                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2080               }
2081               // misidentified protons, kaons
2082               if(pdgProngID[iprong]==2212) {
2083                fillthis="hbpnonIDTot";
2084                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2085               }
2086               if(pdgProngID[iprong]==321) {
2087                fillthis="hbKnonIDTot";
2088                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2089               }
2090             }
2091             break;
2092         default:
2093           break;
2094         }
2095         if(part->PtProng(iprong)>ptmax)ptmax=part->PtProng(iprong);
2096
2097       } //end loop on prongs
2098         //now histograms where we don't need to check identity
2099         fillthis="hbLcpt";
2100         ((TH1F*)listout->FindObject(fillthis))->Fill(part->Pt());
2101         fillthis = "hbDist12toPrim";
2102         ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDist12toPrim());
2103         ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDist23toPrim());
2104         fillthis = "hbSigmaVert";
2105         ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetSigmaVert());
2106         fillthis = "hbDCAs";
2107         Double_t dcas[3];
2108         part->GetDCAs(dcas);
2109         for (Int_t idca=0;idca<3;idca++) ((TH1F*)listout->FindObject(fillthis))->Fill(dcas[idca]);
2110         fillthis = "hbCosPointingAngle";
2111         ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle());
2112         fillthis = "hbDecayLength";
2113         ((TH1F*)listout->FindObject(fillthis))->Fill(part->DecayLength());
2114       Double_t sum2=part->Getd0Prong(0)*part->Getd0Prong(0)+part->Getd0Prong(1)*part->Getd0Prong(1)+part->Getd0Prong(2)*part->Getd0Prong(2);
2115       fillthis = "hbSum2";
2116       ((TH1F*)listout->FindObject(fillthis))->Fill(sum2);
2117       fillthis = "hbptmax";
2118       ((TH1F*)listout->FindObject(fillthis))->Fill(ptmax);
2119
2120      }
2121
2122     } //end mass cut
2123
2124   } // end background case
2125   return;
2126 }
2127