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   if(vtx1==0x0) return;
 821
 822
 823   // load MC particles
 824   if(fReadMC){
 825
 826     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 827     if(!arrayMC) {
 828       printf("AliAnalysisTaskSELambdac::UserExec: MC particles branch not found!\n");
 829       return;
 830     }
 831
 832   // load MC header
 833     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 834     if(!mcHeader) {
 835     printf("AliAnalysisTaskSELambdac::UserExec: MC header branch not found!\n");
 836     return;
 837     }
 838   }
 839
 840   Int_t n3Prong = array3Prong->GetEntriesFast();
 841
 842
 843   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 844     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 845
 846
 847     Bool_t unsetvtx=kFALSE;
 848     if(!d->GetOwnPrimaryVtx()){
 849       d->SetOwnPrimaryVtx(vtx1);
 850       unsetvtx=kTRUE;
 851     }
 852
 853 //    if(d->GetSelectionMap()) {if(!d->HasSelectionBit(AliRDHFCuts::kLcCuts)) continue;}
 854 //   if(d->GetSelectionMap()) if(!d->HasSelectionBit(AliRDHFCuts::kLcPID)) continue;
 855
 856     Int_t isSelectedTracks = fRDCutsProduction->IsSelected(d,AliRDHFCuts::kTracks,aod);
 857     if(!isSelectedTracks) continue;
 858
 859     if (fRDCutsProduction->IsInFiducialAcceptance(d->Pt(),d->Y(4122))) fNentries->Fill(6);
 860
 861     Int_t ptbin=fRDCutsProduction->PtBin(d->Pt());
 862     if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
 863
 864     FillMassHists(aod,d,arrayMC,fRDCutsProduction);
 865     if(fFillVarHists) FillVarHists(d,arrayMC,fRDCutsProduction,fOutputMC,aod);
 866
 867       /*
 868       //start OS analysis
 869       if(labDp<0)fHistOSbkg->Fill(d->InvMassDplus());
 870       fHistOS->Fill(d->InvMassDplus());
 871       */
 872
 873     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 874   }
 875
 876   PostData(1,fOutput);
 877   if(fFillVarHists) PostData(3,fOutputMC);
 878   PostData(4,fNentries);
 879
 880   return;
 881 }
 882
 883
 884
 885 //________________________________________________________________________
 886 void AliAnalysisTaskSELambdac::Terminate(Option_t */*option*/)
 887 {
 888   // Terminate analysis
 889   //
 890   if(fDebug > 1) printf("AnalysisTaskSELambdac: Terminate() \n");
 891
 892   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 893   if (!fOutput) {
 894     printf("ERROR: fOutput not available\n");
 895     return;
 896   }
 897  //fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
 898
 899  TString hisname;
 900   if(fFillNtuple){
 901     fNtupleLambdac = dynamic_cast<TNtuple*>(GetOutputData(3));
 902   }
 903
 904
 905  return;
 906 }
 907
 908 //________________________________________________________________________
 909 Int_t AliAnalysisTaskSELambdac::MatchToMCLambdac(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC) const{
 910   // check if the candidate is a Lambdac decaying in pKpi or in the resonant channels
 911   Int_t lambdacLab[3]={0,0,0};
 912   Int_t pdgs[3]={0,0,0};
 913   for(Int_t i=0;i<3;i++){
 914    AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
 915    Int_t lab=daugh->GetLabel();
 916    if(lab<0) return 0;
 917    AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab);
 918    if(!part) continue;
 919    pdgs[i]=part->GetPdgCode();
 920    Int_t partPdgcode = TMath::Abs(part->GetPdgCode());
 921    if(partPdgcode==211 || partPdgcode==321 || partPdgcode==2212){
 922         Int_t motherLabel=part->GetMother();
 923              if(motherLabel<0) return 0;
 924    AliAODMCParticle *motherPart = (AliAODMCParticle*)arrayMC->At(motherLabel);
 925    if(!motherPart) continue;
 926         Int_t motherPdg = TMath::Abs(motherPart->GetPdgCode());
 927    if(motherPdg==4122) {
 928         if(GetLambdacDaugh(motherPart,arrayMC)){lambdacLab[i]=motherLabel;continue;}
 929               }
 930    if(motherPdg==313 || motherPdg==2224 || motherPdg==3124){
 931          Int_t GmotherLabel=motherPart->GetMother();
 932                if(GmotherLabel<0) return 0;
 933          AliAODMCParticle *GmotherPart = (AliAODMCParticle*)arrayMC->At(GmotherLabel);
 934          if(!GmotherPart) continue;
 935                Int_t GmotherPdg = TMath::Abs(GmotherPart->GetPdgCode());
 936          if(GmotherPdg==4122) {
 937                 if(GetLambdacDaugh(GmotherPart,arrayMC)) {lambdacLab[i]=GmotherLabel;continue;}
 938               }
 939         }
 940      }
 941   }
 942
 943  if(lambdacLab[0]==lambdacLab[1] && lambdacLab[1]==lambdacLab[2]) {return lambdacLab[0];}
 944  return 0;
 945
 946 }
 947 //------------------------
 948 Bool_t AliAnalysisTaskSELambdac::GetLambdacDaugh(AliAODMCParticle *part,TClonesArray *arrayMC) const{
 949  // check if the particle is a lambdac and if its decay mode is the correct one
 950  Int_t numberOfLambdac=0;
 951  if(TMath::Abs(part->GetPdgCode())!=4122) return kFALSE;
 952  Int_t daugh_tmp[2];
 953  daugh_tmp[0]=part->GetDaughter(0);
 954  daugh_tmp[1]=part->GetDaughter(1);
 955  Int_t nDaugh = (Int_t)part->GetNDaughters();
 956  if(nDaugh<2) return kFALSE;
 957  if(nDaugh>3) return kFALSE;
 958  AliAODMCParticle* pdaugh1 = (AliAODMCParticle*)arrayMC->At(part->GetDaughter(0));
 959    if(!pdaugh1) {return kFALSE;}
 960  Int_t number1 = TMath::Abs(pdaugh1->GetPdgCode());
 961  AliAODMCParticle* pdaugh2 = (AliAODMCParticle*)arrayMC->At(part->GetDaughter(1));
 962    if(!pdaugh2) {return kFALSE;}
 963   Int_t number2 = TMath::Abs(pdaugh2->GetPdgCode());
 964
 965   if(nDaugh==3){
 966    Int_t thirdDaugh=part->GetDaughter(1)-1;
 967             AliAODMCParticle* pdaugh3 = (AliAODMCParticle*)arrayMC->At(thirdDaugh);
 968   Int_t number3 = TMath::Abs(pdaugh3->GetPdgCode());
 969    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)) {
 970    numberOfLambdac++;
 971    }
 972   }
 973
 974  if(nDaugh==2){
 975
 976   //Lambda resonant
 977
 978   //Lambda -> p K*0
 979   //
 980   Int_t nfiglieK=0;
 981
 982    if((number1==2212 && number2==313)){
 983      nfiglieK=pdaugh2->GetNDaughters();
 984      if(nfiglieK!=2) return kFALSE;
 985      AliAODMCParticle* pdaughK1 = (AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
 986      AliAODMCParticle* pdaughK2 = (AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
 987      if(!pdaughK1) return kFALSE;
 988      if(!pdaughK2) return kFALSE;
 989      if((TMath::Abs(pdaughK1->GetPdgCode())==211 && TMath::Abs(pdaughK2->GetPdgCode())==321) || (TMath::Abs(pdaughK1->GetPdgCode())==321 && TMath::Abs(pdaughK2->GetPdgCode())==211)) numberOfLambdac++;
 990     }
 991
 992    if((number1==313 && number2==2212)){
 993     nfiglieK=pdaugh1->GetNDaughters();
 994     if(nfiglieK!=2) return kFALSE;
 995     AliAODMCParticle* pdaughK1 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
 996     AliAODMCParticle* pdaughK2 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
 997     if(!pdaughK1) return kFALSE;
 998     if(!pdaughK2) return kFALSE;
 999      if((TMath::Abs(pdaughK1->GetPdgCode())==211 && TMath::Abs(pdaughK2->GetPdgCode())==321) || (TMath::Abs(pdaughK1->GetPdgCode())==321 && TMath::Abs(pdaughK2->GetPdgCode())==211)) numberOfLambdac++;
1000    }
1001
1002    //Lambda -> Delta++ k
1003    Int_t nfiglieDelta=0;
1004    if(number1==321 && number2==2224){
1005     nfiglieDelta=pdaugh2->GetNDaughters();
1006     if(nfiglieDelta!=2) return kFALSE;
1007     AliAODMCParticle *pdaughD1=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
1008     AliAODMCParticle *pdaughD2=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
1009     if(!pdaughD1) return kFALSE;
1010     if(!pdaughD2) return kFALSE;
1011     if((TMath::Abs(pdaughD1->GetPdgCode())==211 && TMath::Abs(pdaughD2->GetPdgCode())==2212) || (TMath::Abs(pdaughD1->GetPdgCode())==2212 && TMath::Abs(pdaughD2->GetPdgCode())==211)) numberOfLambdac++;
1012    }
1013    if(number1==2224 && number2==321){
1014     nfiglieDelta=pdaugh1->GetNDaughters();
1015     if(nfiglieDelta!=2) return kFALSE;
1016     AliAODMCParticle* pdaughD1 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
1017     AliAODMCParticle* pdaughD2 = (AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
1018     if(!pdaughD1) return kFALSE;
1019     if(!pdaughD2) return kFALSE;
1020     if((TMath::Abs(pdaughD1->GetPdgCode())==211 && TMath::Abs(pdaughD2->GetPdgCode())==2212) || (TMath::Abs(pdaughD1->GetPdgCode())==2212 && TMath::Abs(pdaughD2->GetPdgCode())==211)) numberOfLambdac++;
1021    }
1022
1023
1024    //Lambdac -> Lambda(1520) pi
1025    Int_t nfiglieLa=0;
1026    if(number1==3124 && number2==211){
1027     nfiglieLa=pdaugh1->GetNDaughters();
1028     if(nfiglieLa!=2) return kFALSE;
1029     AliAODMCParticle *pdaughL1=(AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(0));
1030     AliAODMCParticle *pdaughL2=(AliAODMCParticle*)arrayMC->At(pdaugh1->GetDaughter(1));
1031     if(!pdaughL1) return kFALSE;
1032     if(!pdaughL2) return kFALSE;
1033     if((TMath::Abs(pdaughL1->GetPdgCode())==321 && TMath::Abs(pdaughL2->GetPdgCode())==2212) || (TMath::Abs(pdaughL1->GetPdgCode())==2212 && TMath::Abs(pdaughL2->GetPdgCode())==321)) numberOfLambdac++;
1034    }
1035    if(number1==211 && number2==3124){
1036     nfiglieLa=pdaugh2->GetNDaughters();
1037     if(nfiglieLa!=2) return kFALSE;
1038     AliAODMCParticle *pdaughL1=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(0));
1039     AliAODMCParticle *pdaughL2=(AliAODMCParticle*)arrayMC->At(pdaugh2->GetDaughter(1));
1040     if(!pdaughL1) return kFALSE;
1041     if(!pdaughL2) return kFALSE;
1042     if((TMath::Abs(pdaughL1->GetPdgCode())==321 && TMath::Abs(pdaughL2->GetPdgCode())==2212) || (TMath::Abs(pdaughL1->GetPdgCode())==2212 && TMath::Abs(pdaughL2->GetPdgCode())==321)) numberOfLambdac++;
1043
1044    }
1045  }
1046
1047  if(numberOfLambdac>0) {return kTRUE;}
1048          return kFALSE;
1049 }
1050 //-----------------------------
1051 Bool_t AliAnalysisTaskSELambdac::IspKpiMC(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC) const{
1052   // Apply MC PID
1053    Int_t lab[3]={0,0,0},pdgs[3]={0,0,0};
1054    for(Int_t i=0;i<3;i++){
1055     AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
1056     lab[i]=daugh->GetLabel();
1057     if(lab[i]<0) return kFALSE;
1058     AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab[i]);
1059     if(!part) return kFALSE;
1060     pdgs[i]=TMath::Abs(part->GetPdgCode());
1061    }
1062
1063    if(pdgs[0]==2212 && pdgs[1]==321 && pdgs[2]==211) return kTRUE;
1064
1065    return kFALSE;
1066 }
1067 //-----------------------------
1068 Bool_t AliAnalysisTaskSELambdac::IspiKpMC(AliAODRecoDecayHF3Prong *d,TClonesArray *arrayMC) const{
1069
1070   // Apply MC PID
1071    Int_t lab[3]={0,0,0},pdgs[3]={0,0,0};
1072    for(Int_t i=0;i<3;i++){
1073     AliAODTrack *daugh=(AliAODTrack*)d->GetDaughter(i);
1074     lab[i]=daugh->GetLabel();
1075     if(lab[i]<0) return kFALSE;
1076     AliAODMCParticle *part= (AliAODMCParticle*)arrayMC->At(lab[i]);
1077     if(!part) return kFALSE;
1078     pdgs[i]=TMath::Abs(part->GetPdgCode());
1079    }
1080
1081    if(pdgs[2]==2212 && pdgs[1]==321 && pdgs[0]==211) {return kTRUE;}
1082
1083    return kFALSE;
1084 }
1085 //--------------------------------------
1086 Bool_t AliAnalysisTaskSELambdac::VertexingKF(AliAODRecoDecayHF3Prong *d,Int_t *pdgs,Double_t field) const{
1087  // apply vertexing KF
1088    Int_t iprongs[3]={0,1,2};
1089    Double_t mass[2]={0.,0.};
1090  //topological constr
1091    AliKFParticle *lambdac=d->ApplyVertexingKF(iprongs,3,pdgs,kTRUE,field,mass);
1092    if(!lambdac) return kFALSE;
1093 //  Double_t probTot=TMath::Prob(Lambdac->GetChi2(),Lambdac->GetNDF());
1094 //  if(probTot<fCutsKF[0]) return kFALSE;
1095   if(lambdac->GetChi2()>fCutsKF[0]) return kFALSE;
1096  //mass constr for K*
1097    Int_t ipRes[2];
1098    Int_t pdgres[2];
1099    mass[0]=0.8961;mass[1]=0.03;
1100    if(TMath::Abs(pdgs[0])==211){
1101     ipRes[0]=0;ipRes[1]=1;
1102     pdgres[0]=pdgs[0];pdgres[1]=321;
1103    }
1104    if(TMath::Abs(pdgs[2])==211){
1105     ipRes[0]=2;ipRes[1]=1;
1106     pdgres[0]=pdgs[2];pdgres[1]=321;
1107    }
1108    AliKFParticle *kappaStar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1109
1110   Double_t probKstar=TMath::Prob(kappaStar->GetChi2(),kappaStar->GetNDF());
1111   if(probKstar>fCutsKF[1]) {
1112     AliAODTrack *esdProng1=(AliAODTrack*)d->GetDaughter(ipRes[0]);
1113     AliAODTrack *esdProng2=(AliAODTrack*)d->GetDaughter(ipRes[1]);
1114     AliKFParticle prong1(*esdProng1,pdgres[0]);
1115     AliKFParticle prong2(*esdProng2,pdgres[1]);
1116     if(kappaStar->GetPt()<fCutsKF[2] && prong1.GetAngle(prong2)>fCutsKF[3]) return kFALSE;
1117   }
1118  //mass constr for Lambda
1119    mass[0]=1.520;mass[1]=0.005;
1120    if(TMath::Abs(pdgs[0])==2212){
1121     ipRes[0]=0;ipRes[1]=1;
1122     pdgres[0]=pdgs[0];pdgres[1]=pdgs[1];
1123    }
1124    if(TMath::Abs(pdgs[2])==2212){
1125     ipRes[0]=2;ipRes[1]=1;
1126     pdgres[0]=pdgs[2];pdgres[1]=pdgs[1];
1127    }
1128    AliKFParticle *lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1129   Double_t probLa=TMath::Prob(lambda1520->GetChi2(),lambda1520->GetNDF());
1130   if(probLa>fCutsKF[4]) {
1131     AliAODTrack *esdProng1=(AliAODTrack*)d->GetDaughter(ipRes[0]);
1132     AliAODTrack *esdProng2=(AliAODTrack*)d->GetDaughter(ipRes[1]);
1133     AliKFParticle prong1(*esdProng1,pdgres[0]);
1134     AliKFParticle prong2(*esdProng2,pdgres[1]);
1135     if(lambda1520->GetPt()<fCutsKF[5] && prong1.GetAngle(prong2)>fCutsKF[6]) return kFALSE;
1136   }
1137  //mass constr for Delta
1138    mass[0]=1.2;mass[1]=0.15;
1139    ipRes[0]=0;ipRes[1]=2;
1140    pdgres[0]=pdgs[0];pdgres[1]=pdgs[2];
1141    AliKFParticle *delta=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1142   Double_t probDelta=TMath::Prob(delta->GetChi2(),delta->GetNDF());
1143   if(probDelta>fCutsKF[7]) {
1144     AliAODTrack *esdProng1=(AliAODTrack*)d->GetDaughter(ipRes[0]);
1145     AliAODTrack *esdProng2=(AliAODTrack*)d->GetDaughter(ipRes[1]);
1146     AliKFParticle prong1(*esdProng1,pdgres[0]);
1147     AliKFParticle prong2(*esdProng2,pdgres[1]);
1148     if(delta->GetPt()<fCutsKF[8] && prong1.GetAngle(prong2)>fCutsKF[9]) return kFALSE;
1149   }
1150  return kTRUE;
1151 }
1152 //-------------------------------------
1153 Int_t AliAnalysisTaskSELambdac::IspiKpResonant(AliAODRecoDecayHF3Prong *d,Double_t field) const{
1154
1155  // apply PID using the resonant channels
1156  //if lambda* -> pk
1157         Double_t mass[2]={1.520,0.005};
1158         Int_t ipRes[2]={1,2};
1159         Int_t pdgres[2]={321,2212};
1160         AliKFParticle *lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1161         Double_t probLa=TMath::Prob(lambda1520->GetChi2(),lambda1520->GetNDF());
1162         if(probLa>0.1) return 1;
1163  //K* -> kpi
1164         mass[0]=0.8961;mass[1]=0.03;
1165         ipRes[0]=0;ipRes[1]=1;
1166         pdgres[0]=211;pdgres[1]=321;
1167         AliKFParticle *kstar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1168         Double_t probKa=TMath::Prob(kstar->GetChi2(),kstar->GetNDF());
1169        if(probKa>0.1) return 2;
1170
1171  return 0;
1172
1173 }
1174 //-------------------------------------
1175 Int_t AliAnalysisTaskSELambdac::IspKpiResonant(AliAODRecoDecayHF3Prong *d,Double_t field) const{
1176
1177  // apply PID using the resonant channels
1178  //if lambda* -> pk
1179         Double_t mass[2]={1.520,0.005};
1180         Int_t ipRes[2]={0,1};
1181         Int_t pdgres[2]={2212,321};
1182         AliKFParticle *lambda1520=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1183         Double_t probLa=TMath::Prob(lambda1520->GetChi2(),lambda1520->GetNDF());
1184         if(probLa>0.1) return 1;
1185  //K* -> kpi
1186         mass[0]=0.8961;mass[1]=0.03;
1187         ipRes[0]=1;ipRes[1]=2;
1188         pdgres[1]=211;pdgres[0]=321;
1189         AliKFParticle *kstar=d->ApplyVertexingKF(ipRes,2,pdgres,kFALSE,field,mass);
1190         Double_t probKa=TMath::Prob(kstar->GetChi2(),kstar->GetNDF());
1191         if(probKa>0.1) return 2;
1192
1193  return 0;
1194
1195 }
1196 //---------------------------
1197 void AliAnalysisTaskSELambdac::FillMassHists(AliAODEvent *aod,AliAODRecoDecayHF3Prong *part, TClonesArray *arrayMC, AliRDHFCutsLctopKpi *cuts){
1198
1199  //if MC PID or no PID, unset PID
1200  if(!fRealPid) cuts->SetUsePID(kFALSE);
1201  Int_t selection=cuts->IsSelected(part,AliRDHFCuts::kCandidate,aod);
1202  if(selection>0){
1203  Int_t iPtBin = -1;
1204  Double_t ptCand = part->Pt();
1205  Int_t index=0;
1206
1207  for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
1208   if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
1209  }
1210
1211
1212  Float_t pdgCode=-2;
1213  Float_t pdgCode1=-1;
1214  Float_t pdgCode2=-1;
1215  Int_t labDp=-1;
1216  Float_t deltaPx=0.;
1217  Float_t deltaPy=0.;
1218  Float_t deltaPz=0.;
1219  Float_t truePt=0.;
1220  Float_t xDecay=0.;
1221  Float_t yDecay=0.;
1222  Float_t zDecay=0.;
1223
1224  if(fReadMC){
1225   labDp = MatchToMCLambdac(part,arrayMC);
1226   if(labDp>0){
1227    AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
1228           AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
1229           AliAODMCParticle *dg1 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(1));
1230           deltaPx=partDp->Px()-part->Px();
1231           deltaPy=partDp->Py()-part->Py();
1232           deltaPz=partDp->Pz()-part->Pz();
1233           truePt=partDp->Pt();
1234           xDecay=dg0->Xv();
1235           yDecay=dg0->Yv();
1236           zDecay=dg0->Zv();
1237           pdgCode=TMath::Abs(partDp->GetPdgCode());
1238           pdgCode1=TMath::Abs(dg0->GetPdgCode());
1239           pdgCode2=TMath::Abs(dg1->GetPdgCode());
1240
1241   }else{
1242    pdgCode=-1;
1243   }
1244  }
1245
1246   Double_t invMasspKpi=-1.;
1247   Double_t invMasspiKp=-1.;
1248   Double_t invMasspKpiLpi=-1.;
1249   Double_t invMasspiKpLpi=-1.;
1250   Double_t invMasspKpiKp=-1.;
1251   Double_t invMasspiKpKp=-1.;
1252   Int_t pdgs[3]={0,0,0};
1253   Double_t field=aod->GetMagneticField();
1254 //apply MC PID
1255   if(fReadMC && fMCPid){
1256
1257   if(IspKpiMC(part,arrayMC)) {
1258    invMasspKpi=part->InvMassLcpKpi();
1259    if(fUseKF){
1260     pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
1261     if(!VertexingKF(part,pdgs,field)) invMasspKpi=-1.;
1262    }
1263   }
1264   if(IspiKpMC(part,arrayMC)) {
1265    invMasspiKp=part->InvMassLcpiKp();
1266    if(fUseKF){
1267     pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
1268     if(!VertexingKF(part,pdgs,field)) invMasspiKp=-1.;
1269    }
1270   }
1271  }
1272
1273  // apply realistic PID
1274   if(fRealPid){
1275    if(selection==1 || selection==3) {
1276     invMasspKpi=part->InvMassLcpKpi();
1277     pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
1278     if(fUseKF){
1279      pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
1280      if(!VertexingKF(part,pdgs,field)) invMasspKpi=-1.;
1281     }
1282    }
1283    if(selection>=2) {
1284     invMasspiKp=part->InvMassLcpiKp();
1285     pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
1286     if(fUseKF){
1287      pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
1288      if(!VertexingKF(part,pdgs,field)) invMasspiKp=-1.;
1289     }
1290    }
1291   }
1292   //apply PID using resonances
1293   if(fResPid && fRealPid){
1294    if((selection==3 || selection==1)) {
1295     if(IspKpiResonant(part,field)==1){
1296      invMasspKpiLpi=part->InvMassLcpKpi();
1297     }
1298     if(IspKpiResonant(part,field)==2){
1299      invMasspKpiKp=part->InvMassLcpKpi();
1300     }
1301    }
1302    if(selection>=2) {
1303     if(IspiKpResonant(part,field)==1){
1304      invMasspiKpLpi=part->InvMassLcpiKp();
1305     }
1306     if(IspiKpResonant(part,field)==2){
1307      invMasspiKpKp=part->InvMassLcpiKp();
1308     }
1309    }
1310   }
1311
1312   // no PID
1313   if(!fResPid && !fRealPid && !fMCPid){
1314    if(selection==2 || selection==3) {
1315     invMasspiKp=part->InvMassLcpiKp();
1316     if(fUseKF){
1317      pdgs[0]=211;pdgs[1]=321;pdgs[2]=2212;
1318      if(!VertexingKF(part,pdgs,field)) invMasspiKp=-1.;
1319     }
1320    }
1321    if(selection==1 || selection==3){
1322     invMasspKpi=part->InvMassLcpKpi();
1323     if(fUseKF){
1324      pdgs[0]=2212;pdgs[1]=321;pdgs[2]=211;
1325      if(!VertexingKF(part,pdgs,field)) invMasspKpi=-1.;
1326      }
1327     }
1328    }
1329
1330    if(invMasspiKp<0. && invMasspKpi<0.) return;
1331
1332    Int_t passTightCuts=0;
1333    if(fAnalysis) passTightCuts=fRDCutsAnalysis->IsSelected(part,AliRDHFCuts::kCandidate,aod);
1334
1335
1336
1337    Float_t tmp[24];
1338       if(fFillNtuple){
1339         tmp[0]=pdgCode;
1340         tmp[1]=deltaPx;
1341         tmp[2]=deltaPy;
1342         tmp[3]=deltaPz;
1343         tmp[4]=truePt;
1344         tmp[5]=xDecay;
1345         tmp[6]=yDecay;
1346         tmp[7]=zDecay;
1347         if(pdgCode1==2212) {tmp[8]=part->PtProng(0);}else{tmp[8]=0.;}
1348         if(pdgCode1==211) {tmp[10]=part->PtProng(0);}else{tmp[10]=0.;}
1349         tmp[9]=part->PtProng(1);
1350         if(pdgCode2==211) {tmp[10]=part->PtProng(2);}else{tmp[10]=0.;}
1351         tmp[11]=part->Pt();
1352         tmp[12]=part->CosPointingAngle();
1353         tmp[13]=part->DecayLength();
1354         tmp[14]=part->Xv();
1355         tmp[15]=part->Yv();
1356         tmp[16]=part->Zv();
1357         if(invMasspiKp>0.) tmp[17]=invMasspiKp;
1358         if(invMasspKpi>0.) tmp[17]=invMasspKpi;
1359         tmp[18]=part->GetSigmaVert();
1360         tmp[19]=part->Getd0Prong(0);
1361         tmp[20]=part->Getd0Prong(1);
1362         tmp[21]=part->Getd0Prong(2);
1363         tmp[22]=part->GetDCA();
1364         tmp[23]=part->Prodd0d0();
1365         fNtupleLambdac->Fill(tmp);
1366         PostData(5,fNtupleLambdac);
1367       }
1368
1369     if(part->Pt()>3.){
1370      if(invMasspiKp>0. && invMasspKpi>0.){
1371       if(invMasspiKp>0.) fhMassPtGreater3->Fill(invMasspiKp,0.5);
1372       if(invMasspKpi>0.) fhMassPtGreater3->Fill(invMasspKpi,0.5);
1373      }else{
1374       if(invMasspiKp>0.) fhMassPtGreater3->Fill(invMasspiKp);
1375       if(invMasspKpi>0.) fhMassPtGreater3->Fill(invMasspKpi);
1376      }
1377      if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1378       if(invMasspiKpLpi>0.) fhMassPtGreater3Lpi->Fill(invMasspiKpLpi,0.5);
1379       if(invMasspKpiLpi>0.) fhMassPtGreater3Lpi->Fill(invMasspKpiLpi,0.5);
1380      }else{
1381       if(invMasspiKpLpi>0.) fhMassPtGreater3Lpi->Fill(invMasspiKpLpi);
1382       if(invMasspKpiLpi>0.) fhMassPtGreater3Lpi->Fill(invMasspKpiLpi);
1383      }
1384      if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1385       if(invMasspiKpKp>0.) fhMassPtGreater3Kp->Fill(invMasspiKpKp,0.5);
1386       if(invMasspKpiKp>0.) fhMassPtGreater3Kp->Fill(invMasspKpiKp,0.5);
1387      }else{
1388       if(invMasspiKpKp>0.) fhMassPtGreater3Kp->Fill(invMasspiKpKp);
1389       if(invMasspKpiKp>0.) fhMassPtGreater3Kp->Fill(invMasspKpiKp);
1390      }
1391      if(passTightCuts>0){
1392       if(invMasspiKp>0. && invMasspKpi>0.){
1393        if(invMasspiKp>0.) fhMassPtGreater3TC->Fill(invMasspiKp,0.5);
1394        if(invMasspKpi>0.) fhMassPtGreater3TC->Fill(invMasspKpi,0.5);
1395       }else{
1396        if(invMasspiKp>0.) fhMassPtGreater3TC->Fill(invMasspiKp);
1397        if(invMasspKpi>0.) fhMassPtGreater3TC->Fill(invMasspKpi);
1398       }
1399      if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1400       if(invMasspiKpLpi>0.) fhMassPtGreater3LpiTC->Fill(invMasspiKpLpi,0.5);
1401       if(invMasspKpiLpi>0.) fhMassPtGreater3LpiTC->Fill(invMasspKpiLpi,0.5);
1402      }else{
1403       if(invMasspiKpLpi>0.) fhMassPtGreater3LpiTC->Fill(invMasspiKpLpi);
1404       if(invMasspKpiLpi>0.) fhMassPtGreater3LpiTC->Fill(invMasspKpiLpi);
1405      }
1406      if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1407       if(invMasspiKpKp>0.) fhMassPtGreater3KpTC->Fill(invMasspiKpKp,0.5);
1408       if(invMasspKpiKp>0.) fhMassPtGreater3KpTC->Fill(invMasspKpiKp,0.5);
1409      }else{
1410       if(invMasspiKpKp>0.) fhMassPtGreater3KpTC->Fill(invMasspiKpKp);
1411       if(invMasspKpiKp>0.) fhMassPtGreater3KpTC->Fill(invMasspKpiKp);
1412      }
1413      }
1414     }
1415    if(part->Pt()>2.){
1416      if(invMasspiKp>0. && invMasspKpi>0.){
1417       if(invMasspiKp>0.) fhMassPtGreater2->Fill(invMasspiKp,0.5);
1418       if(invMasspKpi>0.) fhMassPtGreater2->Fill(invMasspKpi,0.5);
1419      }else{
1420       if(invMasspiKp>0.) fhMassPtGreater2->Fill(invMasspiKp);
1421       if(invMasspKpi>0.) fhMassPtGreater2->Fill(invMasspKpi);
1422      }
1423      if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1424       if(invMasspiKpLpi>0.) fhMassPtGreater2Lpi->Fill(invMasspiKpLpi,0.5);
1425       if(invMasspKpiLpi>0.) fhMassPtGreater2Lpi->Fill(invMasspKpiLpi,0.5);
1426      }else{
1427       if(invMasspiKpLpi>0.) fhMassPtGreater2Lpi->Fill(invMasspiKpLpi);
1428       if(invMasspKpiLpi>0.) fhMassPtGreater2Lpi->Fill(invMasspKpiLpi);
1429      }
1430      if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1431       if(invMasspiKpKp>0.) fhMassPtGreater2Kp->Fill(invMasspiKpKp,0.5);
1432       if(invMasspKpiKp>0.) fhMassPtGreater2Kp->Fill(invMasspKpiKp,0.5);
1433      }else{
1434       if(invMasspiKpKp>0.) fhMassPtGreater2Kp->Fill(invMasspiKpKp);
1435       if(invMasspKpiKp>0.) fhMassPtGreater2Kp->Fill(invMasspKpiKp);
1436      }
1437      if(passTightCuts>0){
1438       if(invMasspiKp>0. && invMasspKpi>0.){
1439        if(invMasspiKp>0.) fhMassPtGreater2TC->Fill(invMasspiKp,0.5);
1440        if(invMasspKpi>0.) fhMassPtGreater2TC->Fill(invMasspKpi,0.5);
1441       }else{
1442        if(invMasspiKp>0.) fhMassPtGreater2TC->Fill(invMasspiKp);
1443        if(invMasspKpi>0.) fhMassPtGreater2TC->Fill(invMasspKpi);
1444       }
1445      if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1446       if(invMasspiKpLpi>0.) fhMassPtGreater2LpiTC->Fill(invMasspiKpLpi,0.5);
1447       if(invMasspKpiLpi>0.) fhMassPtGreater2LpiTC->Fill(invMasspKpiLpi,0.5);
1448      }else{
1449       if(invMasspiKpLpi>0.) fhMassPtGreater2LpiTC->Fill(invMasspiKpLpi);
1450       if(invMasspKpiLpi>0.) fhMassPtGreater2LpiTC->Fill(invMasspKpiLpi);
1451      }
1452      if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1453       if(invMasspiKpKp>0.) fhMassPtGreater2KpTC->Fill(invMasspiKpKp,0.5);
1454       if(invMasspKpiKp>0.) fhMassPtGreater2KpTC->Fill(invMasspKpiKp,0.5);
1455      }else{
1456       if(invMasspiKpKp>0.) fhMassPtGreater2KpTC->Fill(invMasspiKpKp);
1457       if(invMasspKpiKp>0.) fhMassPtGreater2KpTC->Fill(invMasspKpiKp);
1458      }
1459      }
1460     }
1461
1462   if(iPtBin>=0){
1463    index=GetHistoIndex(iPtBin);
1464    if(invMasspiKp>0. && invMasspKpi>0.){
1465     if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp,0.5);
1466     if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi,0.5);
1467    }else{
1468     if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
1469     if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
1470    }
1471    if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1472     if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi,0.5);
1473     if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi,0.5);
1474    }else{
1475     if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi);
1476     if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi);
1477    }
1478    if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1479     if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp,0.5);
1480     if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp,0.5);
1481    }else{
1482     if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp);
1483     if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp);
1484    }
1485    if(passTightCuts>0){
1486     if(invMasspiKp>0. && invMasspKpi>0. && passTightCuts==3){
1487      if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp,0.5);
1488      if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi,0.5);
1489     }else{
1490      if(invMasspiKp>0. && passTightCuts==2) fMassHistTC[index]->Fill(invMasspiKp);
1491      if(invMasspKpi>0. && passTightCuts==1) fMassHistTC[index]->Fill(invMasspKpi);
1492     }
1493    if(invMasspiKpLpi>0. && invMasspKpiLpi>0. && passTightCuts==3){
1494     if(invMasspiKpLpi>0.) fMassHistLpiTC[index]->Fill(invMasspiKpLpi,0.5);
1495     if(invMasspKpiLpi>0.) fMassHistLpiTC[index]->Fill(invMasspKpiLpi,0.5);
1496    }else{
1497     if(invMasspiKpLpi>0. && passTightCuts==2) fMassHistLpiTC[index]->Fill(invMasspiKpLpi);
1498     if(invMasspKpiLpi>0.&& passTightCuts==1) fMassHistLpiTC[index]->Fill(invMasspKpiLpi);
1499    }
1500    if(invMasspiKpKp>0. && invMasspKpiKp>0. && passTightCuts==3){
1501     if(invMasspiKpKp>0.) fMassHistKpTC[index]->Fill(invMasspiKpKp,0.5);
1502     if(invMasspKpiKp>0.) fMassHistKpTC[index]->Fill(invMasspKpiKp,0.5);
1503    }else{
1504     if(invMasspiKpKp>0. && passTightCuts==2) fMassHistKpTC[index]->Fill(invMasspiKpKp);
1505     if(invMasspKpiKp>0.&& passTightCuts==1) fMassHistKpTC[index]->Fill(invMasspKpiKp);
1506    }
1507    }
1508    if(fReadMC){
1509     if(labDp>0) {
1510      index=GetSignalHistoIndex(iPtBin);
1511      if(invMasspiKp>0. && invMasspKpi>0.){
1512       if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp,0.5);
1513       if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi,0.5);
1514      }else{
1515       if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
1516       if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
1517      }
1518     if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1519      if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi,0.5);
1520      if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi,0.5);
1521     }else{
1522      if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi);
1523      if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi);
1524     }
1525     if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1526      if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp,0.5);
1527      if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp,0.5);
1528     }else{
1529      if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp);
1530      if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp);
1531     }
1532      if(passTightCuts>0){
1533       if(invMasspiKp>0. && invMasspKpi>0. && passTightCuts==3){
1534        if(invMasspiKp>0.) fMassHistTC[index]->Fill(invMasspiKp,0.5);
1535        if(invMasspKpi>0.) fMassHistTC[index]->Fill(invMasspKpi,0.5);
1536       }else{
1537        if(invMasspiKp>0. && passTightCuts==2) fMassHistTC[index]->Fill(invMasspiKp);
1538        if(invMasspKpi>0.&& passTightCuts==1) fMassHistTC[index]->Fill(invMasspKpi);
1539       }
1540     if(invMasspiKpLpi>0. && invMasspKpiLpi>0. && passTightCuts==3){
1541      if(invMasspiKpLpi>0.) fMassHistLpiTC[index]->Fill(invMasspiKpLpi,0.5);
1542      if(invMasspKpiLpi>0.) fMassHistLpiTC[index]->Fill(invMasspKpiLpi,0.5);
1543     }else{
1544      if(invMasspiKpLpi>0. && passTightCuts==2) fMassHistLpiTC[index]->Fill(invMasspiKpLpi);
1545      if(invMasspKpiLpi>0.&& passTightCuts==1) fMassHistLpiTC[index]->Fill(invMasspKpiLpi);
1546     }
1547     if(invMasspiKpKp>0. && invMasspKpiKp>0. && passTightCuts==3){
1548      if(invMasspiKpKp>0.) fMassHistKpTC[index]->Fill(invMasspiKpKp,0.5);
1549      if(invMasspKpiKp>0.) fMassHistKpTC[index]->Fill(invMasspKpiKp,0.5);
1550     }else{
1551      if(invMasspiKpKp>0. && passTightCuts==2) fMassHistKpTC[index]->Fill(invMasspiKpKp);
1552      if(invMasspKpiKp>0.&& passTightCuts==1) fMassHistKpTC[index]->Fill(invMasspKpiKp);
1553     }
1554    }
1555
1556    }else{
1557     index=GetBackgroundHistoIndex(iPtBin);
1558     if(invMasspiKp>0. && invMasspKpi>0.){
1559      fMassHist[index]->Fill(invMasspiKp,0.5);
1560      fMassHist[index]->Fill(invMasspKpi,0.5);
1561     }else{
1562      if(invMasspiKp>0.) fMassHist[index]->Fill(invMasspiKp);
1563      if(invMasspKpi>0.) fMassHist[index]->Fill(invMasspKpi);
1564     }
1565     if(invMasspiKpLpi>0. && invMasspKpiLpi>0.){
1566      if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi,0.5);
1567      if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi,0.5);
1568     }else{
1569      if(invMasspiKpLpi>0.) fMassHistLpi[index]->Fill(invMasspiKpLpi);
1570      if(invMasspKpiLpi>0.) fMassHistLpi[index]->Fill(invMasspKpiLpi);
1571     }
1572     if(invMasspiKpKp>0. && invMasspKpiKp>0.){
1573      if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp,0.5);
1574      if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp,0.5);
1575     }else{
1576      if(invMasspiKpKp>0.) fMassHistKp[index]->Fill(invMasspiKpKp);
1577      if(invMasspKpiKp>0.) fMassHistKp[index]->Fill(invMasspKpiKp);
1578     }
1579     if(invMasspiKp>0. && invMasspKpi>0. && passTightCuts==3){
1580      fMassHistTC[index]->Fill(invMasspiKp,0.5);
1581      fMassHistTC[index]->Fill(invMasspKpi,0.5);
1582     }else{
1583      if(invMasspiKp>0. && passTightCuts==2) fMassHistTC[index]->Fill(invMasspiKp);
1584      if(invMasspKpi>0. && passTightCuts==1) fMassHistTC[index]->Fill(invMasspKpi);
1585     }
1586     if(invMasspiKpLpi>0. && invMasspKpiLpi>0. && passTightCuts==3){
1587      if(invMasspiKpLpi>0.) fMassHistLpiTC[index]->Fill(invMasspiKpLpi,0.5);
1588      if(invMasspKpiLpi>0.) fMassHistLpiTC[index]->Fill(invMasspKpiLpi,0.5);
1589     }else{
1590      if(invMasspiKpLpi>0. && passTightCuts==2) fMassHistLpiTC[index]->Fill(invMasspiKpLpi);
1591      if(invMasspKpiLpi>0.&& passTightCuts==1) fMassHistLpiTC[index]->Fill(invMasspKpiLpi);
1592     }
1593     if(invMasspiKpKp>0. && invMasspKpiKp>0. && passTightCuts==3){
1594      if(invMasspiKpKp>0.) fMassHistKpTC[index]->Fill(invMasspiKpKp,0.5);
1595      if(invMasspKpiKp>0.) fMassHistKpTC[index]->Fill(invMasspKpiKp,0.5);
1596     }else{
1597      if(invMasspiKpKp>0. && passTightCuts==2) fMassHistKpTC[index]->Fill(invMasspiKpKp);
1598      if(invMasspKpiKp>0.&& passTightCuts==1) fMassHistKpTC[index]->Fill(invMasspKpiKp);
1599     }
1600    }
1601
1602   }
1603  }
1604 }
1605  return;
1606 }
1607 //-----------------------
1608 void AliAnalysisTaskSELambdac::FillVarHists(AliAODRecoDecayHF3Prong *part, TClonesArray *arrMC, AliRDHFCutsLctopKpi *cuts, TList *listout,AliAODEvent* aod){
1609   //
1610   // function used in UserExec to fill variable histograms analysing MC
1611   //
1612
1613
1614   Int_t pdgDgLctopKpi[3]={2212,321,211};
1615   Int_t lab=-9999;
1616   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)
1617   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)
1618                                                  //3 both (it should never happen...)
1619
1620   if (isSelectedPID==0)fNentries->Fill(7);
1621   if (isSelectedPID==1)fNentries->Fill(8);
1622   if (isSelectedPID==2)fNentries->Fill(9);
1623   if (isSelectedPID==3)fNentries->Fill(10);
1624
1625
1626   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
1627   Double_t minvLcpKpi = part->InvMassLcpKpi();
1628   Double_t minvLcpiKp = part->InvMassLcpiKp();
1629
1630   Double_t invmasscut=0.05;
1631
1632   TString fillthis="";
1633
1634
1635   AliAODTrack *prong0=(AliAODTrack*)part->GetDaughter(0);
1636   AliAODTrack *prong1=(AliAODTrack*)part->GetDaughter(1);
1637   AliAODTrack *prong2=(AliAODTrack*)part->GetDaughter(2);
1638   if (!prong0 || !prong1 || !prong2) {
1639     fNentries->Fill(5);
1640     return;
1641   }
1642
1643   //check pdg of the prongs
1644   Int_t labprong[3];
1645   if(fReadMC){
1646     labprong[0]=prong0->GetLabel();
1647     labprong[1]=prong1->GetLabel();
1648     labprong[2]=prong2->GetLabel();
1649   }
1650   AliAODMCParticle *mcprong=0;
1651
1652   Int_t pdgProng[3]={0,0,0};
1653   Int_t pdgProngID[3]={0,0,0};
1654     if(fReadMC) {
1655      for (Int_t iprong=0;iprong<3;iprong++){
1656       if(labprong[iprong]<0) continue;
1657       mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
1658       pdgProng[iprong]=TMath::Abs(mcprong->GetPdgCode());
1659      }
1660      if(isSelectedPID>0 && fReadMC){
1661       pdgProngID[1]=321;
1662       if(isSelectedPID==1) {pdgProngID[0]=2212;pdgProngID[2]=211;}
1663       if(isSelectedPID==2) {pdgProngID[0]=211;pdgProngID[2]=2212;}
1664      }
1665     } else {
1666       if (isSelectedPID>0) {
1667         pdgProng[1]=321;
1668         if(isSelectedPID==1) {pdgProng[0]=2212;pdgProng[2]=211;}
1669         if(isSelectedPID==2) {pdgProng[0]=211;pdgProng[2]=2212;}
1670       }
1671     }
1672
1673  Int_t selection=cuts->IsSelected(part,AliRDHFCuts::kCandidate,aod);
1674
1675   if((lab>=0 && fReadMC) || (!fReadMC && (isSelectedPID>0)) ){ //signal (from MC or PID)
1676
1677     fillthis="hMass";
1678
1679      if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 4122)
1680         || (!fReadMC && (isSelectedPID>0 && part->Charge()>0))){    //Lc
1681          if ( (pdgProng[0]==2212) && (pdgProng[1]==321) && (pdgProng[2]==211) ) ((TH1F*)listout->FindObject(fillthis))->Fill(minvLcpKpi);
1682          else if ( (pdgProng[0]==211) && (pdgProng[1]==321) && (pdgProng[2]==2212) ) ((TH1F*)listout->FindObject(fillthis))->Fill(minvLcpiKp);
1683     }
1684     else if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == -4122)
1685           || (!fReadMC && (isSelectedPID>0 && part->Charge()<0))){ //anti-Lc
1686       if ( (pdgProng[0]==2212) && (pdgProng[1]==321) && (pdgProng[2]==211) ) ((TH1F*)listout->FindObject(fillthis))->Fill(minvLcpKpi);
1687       else if ( (pdgProng[0]==211) && (pdgProng[1]==321) && (pdgProng[2]==2212) ) ((TH1F*)listout->FindObject(fillthis))->Fill(minvLcpiKp);
1688     }
1689
1690
1691     //apply cut on invariant mass on the pair
1692     if(selection>0){
1693
1694       Double_t ptmax=0.;
1695       for (Int_t iprong=0; iprong<3; iprong++) {
1696         AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(iprong);
1697         if (fReadMC) {
1698          labprong[iprong]=prong->GetLabel();
1699          if(pdgProngID[iprong]==2212) {
1700           fillthis="hpIDTot";
1701           ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1702          }
1703          if(pdgProngID[iprong]==321) {
1704           fillthis="hKIDTot";
1705           ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1706          }
1707          if(pdgProngID[iprong]==211) {
1708           fillthis="hpiIDTot";
1709           ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1710         }
1711        }
1712         AliAODPid *pidObjtrk=prong->GetDetPid();
1713         AliAODPidHF *pidObj=new AliAODPidHF();
1714         Bool_t hasTOF=pidObj->CheckStatus(prong,"TOF");
1715         Bool_t hasTPC=pidObj->CheckStatus(prong,"TPC");
1716         delete pidObj;
1717         Double_t tofSignal=0.;
1718         Double_t dedxTPC=0.;
1719         Double_t momTOF=0.;
1720         Double_t momTPC=0.;
1721         if(hasTOF) {
1722          momTOF = prong->P();
1723          tofSignal=pidObjtrk->GetTOFsignal();
1724         }
1725         if(hasTPC) {
1726          momTPC = pidObjtrk->GetTPCmomentum();
1727          dedxTPC=pidObjtrk->GetTPCsignal();
1728         }
1729         switch (pdgProng[iprong]) {
1730         case 2212:
1731             fillthis="hpTOFSignal";
1732             ((TH1F*)listout->FindObject(fillthis))->Fill(tofSignal);
1733             fillthis="hpTPCSignal";
1734             ((TH1F*)listout->FindObject(fillthis))->Fill(dedxTPC);
1735             fillthis="hpptProng";
1736             ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1737             fillthis="hpd0Prong";
1738             ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
1739             fillthis="hpSignalVspTPC";
1740             ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,dedxTPC);
1741             fillthis="hpSignalVspTOF";
1742             ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,tofSignal);
1743             ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
1744             AliPID::EParticleType typep;
1745             typep=AliPID::EParticleType(4);
1746             if(hasTPC) {
1747              Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typep);
1748              fillthis="hpSigmaVspTPC";
1749              ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,nsigmap);
1750             }
1751             if(hasTOF){
1752              Double_t nsigma=fUtilPid->NumberOfSigmasTOF(prong,typep);
1753              fillthis="hpSigmaVspTOF";
1754              ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,nsigma);
1755             }
1756             if(fReadMC){
1757               // real protons
1758              fillthis="hpRealTot";
1759              ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1760               // id protons
1761               if(pdgProngID[iprong]==2212) {
1762                fillthis="hpIDGood";
1763                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1764               }
1765               // misidentified kaons, pions
1766               if(pdgProngID[iprong]==321) {
1767                fillthis="hKnonIDTot";
1768                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1769               }
1770               if(pdgProngID[iprong]==211) {
1771                fillthis="hpinonIDTot";
1772                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1773               }
1774             }
1775
1776            break;
1777         case 321:
1778             fillthis="hKTOFSignal";
1779             ((TH1F*)listout->FindObject(fillthis))->Fill(tofSignal);
1780             fillthis="hKTPCSignal";
1781             ((TH1F*)listout->FindObject(fillthis))->Fill(dedxTPC);
1782             fillthis="hKptProng";
1783             ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1784             fillthis="hKd0Prong";
1785             ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
1786             fillthis="hKSignalVspTPC";
1787             ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,dedxTPC);
1788             fillthis="hKSignalVspTOF";
1789             ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,tofSignal);
1790             AliPID::EParticleType typek;
1791             typek=AliPID::EParticleType(3);
1792             if(hasTPC) {
1793              Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typek);
1794              fillthis="hKSigmaVspTPC";
1795              ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,nsigmap);
1796             }
1797             if(hasTOF){
1798              Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typek);
1799              fillthis="hKSigmaVspTOF";
1800              ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,nsigma);
1801             }
1802
1803             if(fReadMC){
1804               // real kaons
1805              fillthis="hKRealTot";
1806              ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1807               // id kaons
1808               if(pdgProngID[iprong]==321) {
1809                fillthis="hKIDGood";
1810                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1811               }
1812               // misidentified protons, pions
1813               if(pdgProngID[iprong]==2212) {
1814                fillthis="hpnonIDTot";
1815                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1816               }
1817               if(pdgProngID[iprong]==211) {
1818                fillthis="hpinonIDTot";
1819                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1820               }
1821             }
1822             break;
1823         case 211:
1824             fillthis="hpiTOFSignal";
1825             ((TH1F*)listout->FindObject(fillthis))->Fill(tofSignal);
1826             fillthis="hpiTPCSignal";
1827             ((TH1F*)listout->FindObject(fillthis))->Fill(dedxTPC);
1828             fillthis="hpiptProng";
1829             ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1830             fillthis="hpid0Prong";
1831             ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
1832             fillthis="hpiSignalVspTPC";
1833             ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,dedxTPC);
1834             fillthis="hpiSignalVspTOF";
1835             ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,tofSignal);
1836             AliPID::EParticleType typepi;
1837             typepi=AliPID::EParticleType(2);
1838             if(hasTPC) {
1839              Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typepi);
1840             fillthis="hpiSigmaVspTPC";
1841             ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,nsigmap);
1842             }
1843             if(hasTOF){
1844              Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typepi);
1845              fillthis="hpiSigmaVspTOF";
1846              ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,nsigma);
1847             }
1848
1849             if(fReadMC){
1850               // real pions
1851              fillthis="hpiRealTot";
1852              ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1853               // id pions
1854               if(pdgProngID[iprong]==211) {
1855                fillthis="hpiIDGood";
1856                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1857               }
1858               // misidentified protons, kaons
1859               if(pdgProngID[iprong]==2212) {
1860                fillthis="hpnonIDTot";
1861                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1862               }
1863               if(pdgProngID[iprong]==321) {
1864                fillthis="hKnonIDTot";
1865                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1866               }
1867             }
1868             break;
1869
1870         default:
1871           break;
1872         }
1873         if(part->PtProng(iprong)>ptmax)ptmax=part->PtProng(iprong);
1874
1875       } //end loop on prongs
1876         //now histograms where we don't need to check identity
1877       fillthis = "hDist12toPrim";
1878       ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDist12toPrim());
1879       ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDist23toPrim());
1880       fillthis = "hSigmaVert";
1881       ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetSigmaVert());
1882       fillthis = "hDCAs";
1883       Double_t dcas[3];
1884       part->GetDCAs(dcas);
1885       for (Int_t idca=0;idca<3;idca++) ((TH1F*)listout->FindObject(fillthis))->Fill(dcas[idca]);
1886       fillthis = "hCosPointingAngle";
1887       ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle());
1888       fillthis = "hDecayLength";
1889       ((TH1F*)listout->FindObject(fillthis))->Fill(part->DecayLength());
1890       Double_t sum2=part->Getd0Prong(0)*part->Getd0Prong(0)+part->Getd0Prong(1)*part->Getd0Prong(1)+part->Getd0Prong(2)*part->Getd0Prong(2);
1891       fillthis = "hSum2";
1892       ((TH1F*)listout->FindObject(fillthis))->Fill(sum2);
1893       fillthis = "hptmax";
1894       ((TH1F*)listout->FindObject(fillthis))->Fill(ptmax);
1895       fillthis="hLcpt";
1896       ((TH1F*)listout->FindObject(fillthis))->Fill(part->Pt());
1897
1898     } //end mass cut
1899
1900
1901   } else if( (lab<0) && fReadMC) { // background **ONLY MC**
1902
1903     fillthis="hbMass";
1904
1905      if (part->Charge()>0){    //Lc background
1906          if ( (pdgProng[0]==2212) && (pdgProng[1]==321) && (pdgProng[2]==211) ) ((TH1F*)listout->FindObject(fillthis))->Fill(minvLcpKpi);
1907          else if ( (pdgProng[0]==211) && (pdgProng[1]==321) && (pdgProng[2]==2212) ) ((TH1F*)listout->FindObject(fillthis))->Fill(minvLcpiKp);
1908      }
1909      else if (part->Charge()<0){ //anti-Lc background
1910       if ( (pdgProng[0]==2212) && (pdgProng[1]==321) && (pdgProng[2]==211) ) ((TH1F*)listout->FindObject(fillthis))->Fill(minvLcpKpi);
1911       else if ( (pdgProng[0]==211) && (pdgProng[1]==321) && (pdgProng[2]==2212) ) ((TH1F*)listout->FindObject(fillthis))->Fill(minvLcpiKp);
1912     }
1913
1914
1915     //apply cut on invariant mass on the pair
1916     if(TMath::Abs(minvLcpKpi-mPDG)<invmasscut || TMath::Abs(minvLcpiKp-mPDG)<invmasscut){
1917     if(selection>0){
1918
1919
1920       Double_t ptmax=0.;
1921
1922       for (Int_t iprong=0; iprong<3; iprong++) {
1923         AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(iprong);
1924         if (fReadMC) {
1925         labprong[iprong]=prong->GetLabel();
1926          if(pdgProngID[iprong]==2212) {
1927           fillthis="hbpIDTot";
1928           ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1929          }
1930          if(pdgProngID[iprong]==321) {
1931           fillthis="hbKIDTot";
1932           ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1933          }
1934          if(pdgProngID[iprong]==211) {
1935           fillthis="hbpiIDTot";
1936           ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1937          }
1938         }
1939         AliAODPid *pidObjtrk=prong->GetDetPid();
1940         AliAODPidHF *pidObj=new AliAODPidHF();
1941         Bool_t hasTOF=pidObj->CheckStatus(prong,"TOF");
1942         Bool_t hasTPC=pidObj->CheckStatus(prong,"TPC");
1943         delete pidObj;
1944         Double_t tofSignal=0.;
1945         Double_t dedxTPC=0.;
1946         Double_t momTOF=0.;
1947         Double_t momTPC=0.;
1948         if(hasTOF) {
1949          momTOF = prong->P();
1950          tofSignal=pidObjtrk->GetTOFsignal();
1951         }
1952         if(hasTPC) {
1953          momTPC = pidObjtrk->GetTPCmomentum();
1954          dedxTPC=pidObjtrk->GetTPCsignal();
1955         }
1956
1957         switch (pdgProng[iprong]) {
1958         case 2212:
1959             fillthis="hbpTOFSignal";
1960             ((TH1F*)listout->FindObject(fillthis))->Fill(tofSignal);
1961             fillthis="hbpTPCSignal";
1962             ((TH1F*)listout->FindObject(fillthis))->Fill(dedxTPC);
1963             fillthis="hbpptProng";
1964             ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1965             fillthis="hbpd0Prong";
1966             ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
1967             fillthis="hbpSignalVspTPC";
1968             ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,dedxTPC);
1969             fillthis="hbpSignalVspTOF";
1970             ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,tofSignal);
1971             AliPID::EParticleType typep;
1972             typep=AliPID::EParticleType(4);
1973             if(hasTPC) {
1974              Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typep);
1975              fillthis="hbpSigmaVspTPC";
1976              ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,nsigmap);
1977             }
1978             if(hasTOF){
1979              Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typep);
1980              fillthis="hbpSigmaVspTOF";
1981              ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,nsigma);
1982             }
1983             if(fReadMC){
1984               // real protons
1985              fillthis="hbpRealTot";
1986              ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1987               // id protons
1988               if(pdgProngID[iprong]==2212) {
1989                fillthis="hbpIDGood";
1990                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1991               }
1992               // misidentified kaons, pions
1993               if(pdgProngID[iprong]==321) {
1994                fillthis="hbKnonIDTot";
1995                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
1996               }
1997               if(pdgProngID[iprong]==211) {
1998                fillthis="hbpinonIDTot";
1999                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2000               }
2001             }
2002             break;
2003         case 321:
2004             fillthis="hbKTOFSignal";
2005             ((TH1F*)listout->FindObject(fillthis))->Fill(tofSignal);
2006             fillthis="hbKTPCSignal";
2007             ((TH1F*)listout->FindObject(fillthis))->Fill(dedxTPC);
2008             fillthis="hbKptProng";
2009             ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2010             fillthis="hbKd0Prong";
2011             ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
2012             fillthis="hbKSignalVspTPC";
2013             ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,dedxTPC);
2014             fillthis="hbKSignalVspTOF";
2015             ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,tofSignal);
2016             AliPID::EParticleType typek;
2017             typek=AliPID::EParticleType(3);
2018             if(hasTPC) {
2019              Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typek);
2020              fillthis="hbKSigmaVspTPC";
2021              ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,nsigmap);
2022             }
2023             if(hasTOF){
2024              Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typek);
2025              fillthis="hbKSigmaVspTOF";
2026              ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,nsigma);
2027             }
2028             if(fReadMC){
2029               // real kaons
2030              fillthis="hbKRealTot";
2031              ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2032               // id kaons
2033               if(pdgProngID[iprong]==321) {
2034                fillthis="hbKIDGood";
2035                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2036               }
2037               // misidentified protons, pions
2038               if(pdgProngID[iprong]==2212) {
2039                fillthis="hbpnonIDTot";
2040                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2041               }
2042               if(pdgProngID[iprong]==211) {
2043                fillthis="hbpinonIDTot";
2044                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2045               }
2046             }
2047             break;
2048         case 211:
2049             fillthis="hbpiTOFSignal";
2050             ((TH1F*)listout->FindObject(fillthis))->Fill(tofSignal);
2051             fillthis="hbpiTPCSignal";
2052             ((TH1F*)listout->FindObject(fillthis))->Fill(dedxTPC);
2053             fillthis="hbpiptProng";
2054             ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2055             fillthis="hbpid0Prong";
2056             ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(iprong));
2057             fillthis="hbpiSignalVspTPC";
2058             ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,dedxTPC);
2059             fillthis="hbpiSignalVspTOF";
2060             ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,tofSignal);
2061             AliPID::EParticleType typepi;
2062             typepi=AliPID::EParticleType(2);
2063             if(hasTPC) {
2064              Double_t nsigmap = fUtilPid->NumberOfSigmasTPC(prong,typepi);
2065              fillthis="hbpiSigmaVspTPC";
2066              ((TH2F*)listout->FindObject(fillthis))->Fill(momTPC,nsigmap);
2067             }
2068             if(hasTOF){
2069              Double_t nsigma = fUtilPid->NumberOfSigmasTOF(prong,typepi);
2070              fillthis="hbpiSigmaVspTOF";
2071              ((TH2F*)listout->FindObject(fillthis))->Fill(momTOF,nsigma);
2072             }
2073             if(fReadMC){
2074               // real pions
2075              fillthis="hbpiRealTot";
2076              ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2077               // id pions
2078               if(pdgProngID[iprong]==211) {
2079                fillthis="hbpiIDGood";
2080                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2081               }
2082               // misidentified protons, kaons
2083               if(pdgProngID[iprong]==2212) {
2084                fillthis="hbpnonIDTot";
2085                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2086               }
2087               if(pdgProngID[iprong]==321) {
2088                fillthis="hbKnonIDTot";
2089                ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(iprong));
2090               }
2091             }
2092             break;
2093         default:
2094           break;
2095         }
2096         if(part->PtProng(iprong)>ptmax)ptmax=part->PtProng(iprong);
2097
2098       } //end loop on prongs
2099         //now histograms where we don't need to check identity
2100         fillthis="hbLcpt";
2101         ((TH1F*)listout->FindObject(fillthis))->Fill(part->Pt());
2102         fillthis = "hbDist12toPrim";
2103         ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDist12toPrim());
2104         ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDist23toPrim());
2105         fillthis = "hbSigmaVert";
2106         ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetSigmaVert());
2107         fillthis = "hbDCAs";
2108         Double_t dcas[3];
2109         part->GetDCAs(dcas);
2110         for (Int_t idca=0;idca<3;idca++) ((TH1F*)listout->FindObject(fillthis))->Fill(dcas[idca]);
2111         fillthis = "hbCosPointingAngle";
2112         ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle());
2113         fillthis = "hbDecayLength";
2114         ((TH1F*)listout->FindObject(fillthis))->Fill(part->DecayLength());
2115       Double_t sum2=part->Getd0Prong(0)*part->Getd0Prong(0)+part->Getd0Prong(1)*part->Getd0Prong(1)+part->Getd0Prong(2)*part->Getd0Prong(2);
2116       fillthis = "hbSum2";
2117       ((TH1F*)listout->FindObject(fillthis))->Fill(sum2);
2118       fillthis = "hbptmax";
2119       ((TH1F*)listout->FindObject(fillthis))->Fill(ptmax);
2120
2121      }
2122
2123     } //end mass cut
2124
2125   } // end background case
2126   return;
2127 }
2128