PWGHF/vertexingHF/upgrade/AliAnalysisTaskCountLcEta.cxx

   1 //#####################################################
   2 //#                                                   #
   3 //#          Analysis Task for Lc analysis on ESD     #
   4 //#Authors: C. Bianchin (Utrecht University)          #
   5 //#         and R. Romita (Univ of Liverpool,         #
   6 //#         Daresbury Lab),                           #
   7 //#         based on a class                          #
   8 //#         by MinJung Kweon, Universitaet Heidelberg #
   9 //#                                                   #
  10 //#####################################################
  11
  12 #include "TChain.h"
  13 #include "TTree.h"
  14 #include "TList.h"
  15 #include "TH2F.h"
  16 #include "TF1.h"
  17 #include "TCanvas.h"
  18 #include "TROOT.h"
  19 #include "TNtuple.h"
  20 #include "TLorentzVector.h"
  21 #include "TParticle.h"
  22 #include "TDatabasePDG.h"
  23 #include <TStopwatch.h>
  24 #include <TClonesArray.h>
  25 #include <exception>
  26
  27 #include "AliAnalysisTaskSE.h"
  28 #include "AliAnalysisManager.h"
  29
  30 #include "AliLog.h"
  31 #include "AliStack.h"
  32 #include "AliESDEvent.h"
  33 #include "AliESDInputHandler.h"
  34 #include "AliAODEvent.h"
  35 #include "AliAODInputHandler.h"
  36 #include "AliMCEvent.h"
  37 #include "AliMCEventHandler.h"
  38
  39 #include "AliAnalysisTaskCountLcEta.h"
  40
  41
  42 ClassImp(AliAnalysisTaskCountLcEta) // adding the class to ROOT
  43
  44
  45 //__________________________________________________________________
  46 AliAnalysisTaskCountLcEta::AliAnalysisTaskCountLcEta(const char *name, Int_t ncuts,Double_t *cuts)
  47 : AliAnalysisTaskSE(name)
  48   , fESD(0)
  49   , fAOD(0)
  50   , fAnalysisType("ESD")
  51   , fEvt(0)
  52   , fOutList(0)
  53   , fEnableMCQA(kTRUE)
  54   , fhNevt(0)
  55   , fEtaAbs(0.9)
  56   , fEtaAbsMax(1.5)
  57   , fFillBkg(0)
  58   , fNcuts(ncuts)
  59   , fCuts(cuts)
  60   , fCutNames(0)
  61   , fLooserPtTrack(0)
  62   , fInvMassCut(0.024)
  63   , fThreesigmas(0)
  64 {
  65   // Standard constructor
  66
  67   // Define input and output slots here
  68   DefineInput(0, TChain::Class()); // Input slot #0 works with a TChain
  69   DefineOutput(1, TList::Class()); // Output slot #0 writes into a TList container for mcQA
  70   const Int_t ncutsrightnow=3;
  71   fCutNames=new TString[ncutsrightnow];
  72   if(fNcuts!=ncutsrightnow) {
  73     Printf("ERROR!! Given %d cuts while %d are expected! Fix the class",fNcuts,ncutsrightnow);
  74     return;
  75   }
  76   fCutNames[0]="ptpi";
  77   fCutNames[1]="ptK";
  78   fCutNames[2]="ptp";
  79   Double_t pt=9999999.;
  80   for (Int_t i=0; i<fNcuts; i++){
  81     if(fCutNames[i].Contains("pt") && pt > fCuts[i]) pt=fCuts[i];
  82   }
  83   fLooserPtTrack=pt;
  84   Printf("INFO: Looser pt cuts = %f",fLooserPtTrack);
  85
  86   Double_t mLc=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
  87   fThreesigmas=new Double_t[2];
  88   fThreesigmas[0]=mLc-fInvMassCut;
  89   fThreesigmas[1]=mLc+fInvMassCut;
  90
  91 }
  92
  93 //__________________________________________________________________
  94 AliAnalysisTaskCountLcEta::AliAnalysisTaskCountLcEta(): AliAnalysisTaskSE()
  95   , fESD(0)
  96   , fAOD(0)
  97   , fAnalysisType("ESD")
  98   , fEvt(0)
  99   , fOutList(0)
 100   , fEnableMCQA(kTRUE)
 101   , fhNevt(0)
 102   , fEtaAbs(0.9)
 103   , fEtaAbsMax(1.5)
 104   , fFillBkg(0)
 105   , fNcuts(3)
 106   , fCuts(0)
 107   , fCutNames(0)
 108   , fLooserPtTrack(0)
 109   , fInvMassCut(0.024)
 110   , fThreesigmas(0)
 111   {
 112      //default constructor
 113      Double_t mLc=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
 114      fThreesigmas=new Double_t[2];
 115      fThreesigmas[0]=mLc-fInvMassCut;
 116      fThreesigmas[1]=mLc+fInvMassCut;
 117
 118   }
 119
 120 //__________________________________________________________________
 121 AliAnalysisTaskCountLcEta::AliAnalysisTaskCountLcEta(const AliAnalysisTaskCountLcEta& source): AliAnalysisTaskSE()
 122   , fESD(source.fESD)
 123   , fAOD(source.fAOD)
 124   , fAnalysisType(source.fAnalysisType)
 125   , fEvt(source.fEvt)
 126   , fOutList(source.fOutList)
 127   , fEnableMCQA(source.fEnableMCQA)
 128   , fhNevt(source.fhNevt)
 129   , fEtaAbs(source.fEtaAbs)
 130   , fEtaAbsMax(source.fEtaAbsMax)
 131   , fFillBkg(source.fFillBkg)
 132   , fNcuts(source.fNcuts)
 133   , fCuts(source.fCuts)
 134   , fCutNames(source.fCutNames)
 135   , fLooserPtTrack(source.fLooserPtTrack)
 136   , fInvMassCut(source.fInvMassCut)
 137   , fThreesigmas(source.fThreesigmas)
 138   {
 139    //copy constructor
 140
 141   }
 142
 143 //__________________________________________________________________
 144 void AliAnalysisTaskCountLcEta::UserCreateOutputObjects() {
 145   // Create histograms. Called once
 146
 147   //printf("CreateOutputObjects of task %s\n", GetName());
 148
 149   fOutList = new TList();
 150   fOutList->SetOwner();
 151   // Open a output file
 152   //OpenFile(0);
 153
 154   fhNevt=new TH1F("fhNevt","Number of events",1,-0.5,0.5);
 155   fOutList->Add(fhNevt);
 156   //create histograms for Lambdac
 157
 158   TString hname="hLc3Prongs";
 159   TH1F *hLc3Prongs=new TH1F(hname,"Pt of generated Lambdac in 3 prongs;p_{T} (GeV/c)",100,0.,20.);
 160   hLc3Prongs->Sumw2();
 161   fOutList->Add(hLc3Prongs);
 162
 163   hname="hPLc3Prongs";
 164   TH1F *hPLc3Prongs=new TH1F(hname,"P of generated Lambdac in 3 prongs;p (GeV/c)",100,0.,20.);
 165   hPLc3Prongs->Sumw2();
 166   fOutList->Add(hPLc3Prongs);
 167
 168   hname="hPtpi";
 169   TH1F *hPtpi=new TH1F(hname,"Pt of Lc pi;p_{T} (GeV/c)",100,0.,20.);
 170   hPtpi->Sumw2();
 171   fOutList->Add(hPtpi);
 172   hname="hPpi";
 173   TH1F *hPpi=new TH1F(hname,"P of Lc pi;p (GeV/c)",100,0.,20.);
 174   hPpi->Sumw2();
 175   fOutList->Add(hPpi);
 176
 177   hname="hPtp";
 178   TH1F *hPtp=new TH1F(hname,"Pt of Lc p;p_{T} (GeV/c)",100,0.,20.);
 179   hPtp->Sumw2();
 180   fOutList->Add(hPtp);
 181   hname="hPp";
 182   TH1F *hPp=new TH1F(hname,"P of Lc p;p (GeV/c)",100,0.,20.);
 183   hPp->Sumw2();
 184   fOutList->Add(hPp);
 185
 186   hname="hPtK";
 187   TH1F *hPtK=new TH1F(hname,"Pt of Lc K;p_{T} (GeV/c)",100,0.,20.);
 188   hPtK->Sumw2();
 189   fOutList->Add(hPtK);
 190   hname="hPK";
 191   TH1F *hPK=new TH1F(hname,"P of Lc K;p (GeV/c)",100,0.,20.);
 192   hPK->Sumw2();
 193   fOutList->Add(hPK);
 194
 195   hname="hEtaLc";
 196   TH1F *hEtaLc=new TH1F(hname,"Eta of Lc",400,-20.,20.);
 197   hEtaLc->Sumw2();
 198   fOutList->Add(hEtaLc);
 199
 200   hname="hYLc";
 201   TH1F* hYLc=new TH1F(hname,"Rapidity (y) of Lc",400,-20.,20.);
 202   hYLc->Sumw2();
 203   fOutList->Add(hYLc);
 204
 205
 206   hname="hLc3ProngsInEta";
 207   TH1F *hLc3ProngsInEta=new TH1F(hname,Form("Pt of generated Lambdac in 3 prongs within %.1f ;p_{T} (GeV/c)",fEtaAbs),100,0.,20.);
 208   hLc3ProngsInEta->Sumw2();
 209   fOutList->Add(hLc3ProngsInEta);
 210
 211   hname="hLc3DaughInEta";
 212   TH1F *hLc3DaughInEta=new TH1F(hname,Form("Pt of generated Lambdac with 3 daughters within %.1f;p_{T} (GeV/c)",fEtaAbs),100,0.,20.);
 213   hLc3DaughInEta->Sumw2();
 214   fOutList->Add(hLc3DaughInEta);
 215
 216   hname="hLc2DaughInEta";
 217   TH1F *hLc2DaughInEta=new TH1F(hname,Form("Pt of generated Lambdac with 2 daughters within %.1f;p_{T} (GeV/c)",fEtaAbs),100,0.,20.);
 218   hLc2DaughInEta->Sumw2();
 219   fOutList->Add(hLc2DaughInEta);
 220
 221   hname="hLc1DaughInEta";
 222   TH1F *hLc1DaughInEta=new TH1F(hname,Form("Pt of generated Lambdac with 1 daughter within %.1f;p_{T} (GeV/c)",fEtaAbs),100,0.,20.);
 223   hLc1DaughInEta->Sumw2();
 224   fOutList->Add(hLc1DaughInEta);
 225
 226   TH1F* hLcpiInEta=new TH1F("hLcpiInEta", Form("Pt of generated Lambdac with #pi%.1f < |#eta| < %.1f;p_{T} (GeV/c)",fEtaAbs,fEtaAbsMax),100,0.,20.);
 227   hLcpiInEta->Sumw2();
 228   fOutList->Add(hLcpiInEta);
 229
 230   TH1F* hLcpInEta=new TH1F("hLcpInEta", Form("Pt of generated Lambdac with p %.1f < |#eta| < %.1f;p_{T} (GeV/c)",fEtaAbs,fEtaAbsMax),100,0.,20.);
 231   hLcpInEta->Sumw2();
 232   fOutList->Add(hLcpInEta);
 233
 234   TH1F* hLcKInEta=new TH1F("hLcKInEta", Form("Pt of generated Lambdac with K %.1f < |#eta| < %.1f;p_{T} (GeV/c)",fEtaAbs,fEtaAbsMax),100,0.,20.);
 235   hLcKInEta->Sumw2();
 236   fOutList->Add(hLcKInEta);
 237
 238   TH1F* hLcpiKInEta=new TH1F("hLcpiKInEta", Form("Pt of generated Lambdac with #pi and K %.1f < |#eta| < %.1f;p_{T} (GeV/c)",fEtaAbs,fEtaAbsMax),100,0.,20.);
 239   hLcpiKInEta->Sumw2();
 240   fOutList->Add(hLcpiKInEta);
 241
 242   TH1F* hLcpipInEta=new TH1F("hLcpipInEta", Form("Pt of generated Lambdac with #pi and p %.1f < |#eta| < %.1f;p_{T} (GeV/c)",fEtaAbs,fEtaAbsMax),100,0.,20.);
 243   hLcpipInEta->Sumw2();
 244   fOutList->Add(hLcpipInEta);
 245
 246   TH1F* hLcKpInEta=new TH1F("hLcKpInEta", Form("Pt of generated Lambdac with K and p %.1f < |#eta| < %.1f;p_{T} (GeV/c)",fEtaAbs,fEtaAbsMax),100,0.,20.);
 247   hLcKpInEta->Sumw2();
 248   fOutList->Add(hLcKpInEta);
 249
 250   TH1F* hLcpKpiInEta=new TH1F("hLcpKpiInEta", Form("Pt of generated Lambdac with #pi K and p %.1f < |#eta| < %.1f;p_{T} (GeV/c)",fEtaAbs,fEtaAbsMax),100,0.,20.);
 251   hLcpKpiInEta->Sumw2();
 252   fOutList->Add(hLcpKpiInEta);
 253
 254   TH1F* hRejection=new TH1F("hRejection","Reason of track rejection",7,-0.5,6.5);
 255   hRejection->GetXaxis()->SetBinLabel(1,"not primary");
 256   hRejection->GetXaxis()->SetBinLabel(2,"out of beam pipe");
 257   hRejection->GetXaxis()->SetBinLabel(3,"p_{T} cut");
 258   hRejection->GetXaxis()->SetBinLabel(4,Form("|#eta|>%.1f",fEtaAbs));
 259   hRejection->GetXaxis()->SetBinLabel(5,"Track Selected!");
 260   hRejection->GetXaxis()->SetBinLabel(6,"Candidate Sel (pt cuts)");
 261   hRejection->GetXaxis()->SetBinLabel(7,"Candidate Sel (also mass)");
 262   hRejection->GetXaxis()->SetNdivisions(1,kFALSE);
 263   fOutList->Add(hRejection);
 264   //background
 265   if(fFillBkg){
 266         TH2F* hPtEtaBkg=new TH2F("hPtEtaBkg","P_{T} distribution of background tracks;p_{T} (GeV/c);#eta",100,0.,20.,40,-10.,10.);
 267         fOutList->Add(hPtEtaBkg);
 268
 269         TH1F* hPtEtaMCCandBMid=new TH1F("hPtEtaMCCandBMid","p_{T} distribution of background candidates |#eta| <0.8;p_{T} (GeV/c)",100,0.,20.);
 270         hPtEtaMCCandBMid->Sumw2();
 271         fOutList->Add(hPtEtaMCCandBMid);
 272
 273         TH1F* hPtEtaMCCandBUpg=new TH1F("hPtEtaMCCandBUpg","p_{T} distribution of background candidates |#eta| <1.5;p_{T} (GeV/c)",100,0.,20.);
 274         hPtEtaMCCandBUpg->Sumw2();
 275         fOutList->Add(hPtEtaMCCandBUpg);
 276
 277         TH1F* hMassEtaMCCandBMid=new TH1F("hMassEtaMCCandBMid","Invariant mass distribution of background candidates |#eta| <0.8;inv mass (GeV)",400,2.261,2.309);
 278         hMassEtaMCCandBMid->Sumw2();
 279         fOutList->Add(hMassEtaMCCandBMid);
 280
 281         TH1F* hMassEtaMCCandBUpg=new TH1F("hMassEtaMCCandBUpg","Invariant mass distribution of background candidates |#eta| <1.5;inv mass (GeV)",400,2.261,2.309);
 282         hMassEtaMCCandBUpg->Sumw2();
 283         fOutList->Add(hMassEtaMCCandBUpg);
 284
 285
 286  }
 287
 288 PostData(1,fOutList);
 289
 290 }
 291
 292
 293 //__________________________________________________________________
 294 void AliAnalysisTaskCountLcEta::UserExec(Option_t *) {
 295   // Main loop. Called for each event
 296   fESD=dynamic_cast<AliESDEvent*> (InputEvent());
 297   if (!fESD) { Printf("ERROR: fESD not available"); return;}
 298   // MC Event
 299   AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 300   if (!mcHandler) { printf("ERROR: Could not retrieve MC event handler"); return;}
 301
 302   AliMCEvent* mcEvent = mcHandler->MCEvent();
 303   if (!mcEvent) { Printf("ERROR: Could not retrieve MC event"); return;}
 304   printf("MC Event %p",mcEvent);
 305   AliStack* stack = mcEvent->Stack();
 306   if (!stack) { printf( "Stack not available"); return;}
 307   const AliVVertex *vtx=mcEvent->GetPrimaryVertex();
 308   Double_t position[3]={0.,0.,0.};
 309   vtx->GetXYZ(position);
 310
 311   // fill to count total number of event
 312   fhNevt->Fill(0);
 313
 314   // run MC QA
 315   if (fEnableMCQA) {
 316
 317      //histograms
 318      TH2F* hPtEtaBkg=(TH2F*)fOutList->FindObject("hPtEtaBkg");
 319
 320     //Int_t nPrims = stack->GetNprimary();
 321     Int_t nMCTracks = stack->GetNtrack();
 322
 323     Printf("Loop on %d tracks, %d combinations max \n",nMCTracks, nMCTracks*(nMCTracks-1)*(nMCTracks-2));
 324     TStopwatch time;
 325     time.Start();
 326     // loop over particle in the stack and save the selected ones in an array
 327     TClonesArray arrPartSelpos("TParticle",nMCTracks);
 328     TClonesArray arrPartSelneg("TParticle",nMCTracks);
 329     Int_t nPartSelpos=0, nPartSelneg=0, nLc=0;
 330     for (Int_t istack = 0; istack < nMCTracks; istack++){
 331       TParticle* particle=(TParticle*)stack->Particle(istack);
 332       if(!particle) continue;
 333       Int_t ch=particle->GetPdgCode();
 334       Int_t pdgcode=TMath::Abs(ch);
 335
 336       if(pdgcode==4122) { //signal
 337          FillHistosL(particle,stack);
 338          nLc++;
 339          continue;
 340       }
 341       //background
 342       Bool_t selected=SelectTrack(particle,kTRUE);
 343       if(!selected) continue;
 344
 345       if(ch>0) {
 346          new(arrPartSelpos[nPartSelpos]) TParticle(*particle);
 347          nPartSelpos++;
 348       }
 349       if(ch<0) {
 350          new(arrPartSelneg[nPartSelneg]) TParticle(*particle);
 351          nPartSelneg++;
 352       }
 353
 354     }
 355     Printf("INFO: %d Lc, %d positive and %d negative background particle selected", nLc, nPartSelpos, nPartSelneg);
 356
 357     // loop over primary particles for quark and heavy hadrons
 358     for (Int_t igen = 0; igen < nPartSelpos; igen++){ //nMCTracks
 359       //TParticle* particle=(TParticle*)stack->Particle(igen);
 360       TParticle* particle=(TParticle*)arrPartSelpos.UncheckedAt(igen);
 361       Int_t ch1=particle->GetPdgCode();
 362       Int_t pdgcode=TMath::Abs(ch1);
 363       if(ch1>0) ch1=+1;
 364       if(ch1<0) ch1=-1;
 365
 366
 367
 368       if(fFillBkg && particle->IsPrimary() && (pdgcode==2212 || pdgcode==321 || pdgcode==211))
 369          hPtEtaBkg->Fill(particle->Pt(),particle->Eta());
 370
 371
 372       Bool_t selected=SelectTrack(particle,kTRUE);
 373       if(!selected) continue;
 374       //Printf("Selected First loop");
 375       for(Int_t j=0;fFillBkg && j<nPartSelneg;j++){//second loop
 376
 377         if(igen==j) continue;
 378         TParticle* part2=(TParticle*)arrPartSelneg.UncheckedAt(j);//(TParticle*)stack->Particle(j);
 379
 380
 381         for(Int_t k=0;fFillBkg && k<nPartSelpos;k++){//third loop on pos
 382
 383           if(igen==k || j==k) continue;
 384
 385           TParticle* part3=(TParticle*)arrPartSelpos.UncheckedAt(k);//(TParticle*)stack->Particle(k);
 386
 387           Float_t eta1=particle->Eta();
 388           Float_t eta2=part2->Eta();
 389           Float_t eta3=part3->Eta();
 390           Float_t etamax=TMath::Abs(eta1); if(TMath::Abs(eta2)>etamax) etamax=TMath::Abs(eta2); if(TMath::Abs(eta3)>etamax) etamax=TMath::Abs(eta3);
 391           if(etamax>fEtaAbsMax) continue;
 392           FillHistogramsBackgroundCandidates(particle, part2, part3, etamax);
 393
 394         }//end third loop on pos
 395
 396         for(Int_t k=0;fFillBkg && k<nPartSelneg;k++){//third loop on neg
 397
 398           if(igen==k || j==k) continue;
 399
 400           TParticle* part3=(TParticle*)arrPartSelneg.UncheckedAt(k);//(TParticle*)stack->Particle(k);
 401
 402           Float_t eta1=particle->Eta();
 403           Float_t eta2=part2->Eta();
 404           Float_t eta3=part3->Eta();
 405           Float_t etamax=TMath::Abs(eta1); if(TMath::Abs(eta2)>etamax) etamax=TMath::Abs(eta2); if(TMath::Abs(eta3)>etamax) etamax=TMath::Abs(eta3);
 406           if(etamax>fEtaAbsMax) continue;
 407           FillHistogramsBackgroundCandidates(particle, part2, part3, etamax);
 408
 409         }//end third loop on neg
 410
 411       }//end second loop
 412
 413     }//end loop on generated tracks
 414     time.Stop();
 415     time.Print();
 416     time.Start();
 417
 418     arrPartSelpos.Delete();
 419     arrPartSelneg.Delete();
 420
 421
 422   } // end of MC QA loop
 423
 424
 425   fEvt++;               // event number
 426
 427   PostData(1, fOutList);
 428 }
 429
 430 Bool_t AliAnalysisTaskCountLcEta::SelectTrack(TParticle *p,Bool_t fillh){
 431
 432   TH1F* hrej=(TH1F*)fOutList->FindObject("hRejection");
 433
 434   if(!p->IsPrimary()) {
 435     if(fillh) hrej->Fill(0);
 436     return kFALSE;
 437   }
 438   Float_t bpipe=2.8; //cm
 439   Float_t zvtxdist=1.; //cm
 440
 441
 442   Double_t vx=p->Vx(),vy=p->Vy(),vz=p->Vz();
 443   //Printf("Production point %f, %f", vx,vy);
 444   if((vx*vx + vy*vy) > bpipe*bpipe || vz > zvtxdist) {
 445     if(fillh) hrej->Fill(1);
 446     return kFALSE;
 447   }
 448   Double_t pt=p->Pt();
 449   if(pt<fLooserPtTrack){
 450     if(fillh) hrej->Fill(2);
 451     return kFALSE;
 452   }
 453   Double_t eta=TMath::Abs(p->Eta());
 454   if(eta>fEtaAbs) {
 455     if(fillh) hrej->Fill(3);
 456     if(eta>fEtaAbsMax) return kFALSE;
 457   }
 458   hrej->Fill(4);
 459   return kTRUE;
 460 }
 461
 462 Bool_t AliAnalysisTaskCountLcEta::SelectTracksForCandidate(TParticle* pion, TParticle* kaon, TParticle* proton){
 463
 464   if(pion->Pt()<fCuts[0]) return kFALSE;
 465   if(kaon->Pt()<fCuts[1]) return kFALSE;
 466   if(proton->Pt()<fCuts[2]) return kFALSE;
 467    TH1F* hrej=(TH1F*)fOutList->FindObject("hRejection");
 468   hrej->Fill(5);
 469
 470   return kTRUE;
 471 }
 472
 473 Double_t AliAnalysisTaskCountLcEta::InvMass(TParticle *p1p,TParticle *pn,TParticle *p2p,
 474         TLorentzVector *&candp1ppnp2p){
 475
 476   // TStopwatch watch;
 477   // watch.Start();
 478   //the TLorentzVector is created with NEW, remember to delete it!!
 479   Double_t pxp1=p1p->Px(),pyp1=p1p->Py(),pzp1=p1p->Pz();
 480   Double_t pxp2=p2p->Px(),pyp2=p2p->Py(),pzp2=p2p->Pz();
 481   Double_t pxpn=pn->Px(),pypn=pn->Py(),pzpn=pn->Pz();
 482   //Printf("pxp1 = %f, pyp1 = %f,pzp1=%f ",pxp1,pyp1,pzp1);
 483   Double_t p[3];
 484   p[0]=pxp1+pxpn+pxp2;
 485   p[1]=pyp1+pypn+pyp2;
 486   p[2]=pzp1+pzpn+pzp2;
 487   Double_t masspi=TDatabasePDG::Instance()->GetParticle(211)->Mass();
 488   Double_t massp= TDatabasePDG::Instance()->GetParticle(2212)->Mass();
 489   Double_t massK= TDatabasePDG::Instance()->GetParticle(321)->Mass();
 490   Double_t energyP1=TMath::Sqrt(massp*massp+(p1p->P()*p1p->P()));
 491   Double_t energyPn=TMath::Sqrt(massK*massK+(pn->P()*pn->P()));
 492   Double_t energyP2=TMath::Sqrt(masspi*masspi+(p2p->P()*p2p->P()));
 493
 494
 495   Double_t energy=energyP1 + energyPn +energyP2;
 496   Double_t p2=p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
 497
 498   Double_t invmass=TMath::Sqrt(energy*energy-p2);
 499   candp1ppnp2p=new TLorentzVector(p[0],p[1],p[2],energy);
 500   // watch.Stop();
 501   // watch.Print();
 502   return invmass;
 503
 504 }
 505
 506 void AliAnalysisTaskCountLcEta::FillHistogramsBackgroundCandidates(TParticle *p1,TParticle *p2,TParticle *p3, Double_t etamax){
 507
 508     //histograms
 509     TH1F* hPtEtaCandMid=(TH1F*)fOutList->FindObject("hPtEtaMCCandBMid");
 510     TH1F* hMassEtaCandMid=(TH1F*)fOutList->FindObject("hMassEtaMCCandBMid");
 511     TH1F* hPtEtaCandUpg=(TH1F*)fOutList->FindObject("hPtEtaMCCandBUpg");
 512     TH1F* hMassEtaCandUpg=(TH1F*)fOutList->FindObject("hMassEtaMCCandBUpg");
 513     TH1F* hrej=(TH1F*)fOutList->FindObject("hRejection");
 514
 515    //candidates
 516    Bool_t hyppKpi=kTRUE;
 517    hyppKpi=SelectTracksForCandidate(p1,p2,p3);
 518    Bool_t hyppiKp=kTRUE;
 519    hyppiKp=SelectTracksForCandidate(p3,p2,p1);
 520    if(!hyppKpi && !hyppiKp) return;
 521    Double_t invmasspKpi=0,invmasspiKp=0;
 522    //Printf("Selected Candidates");
 523    TLorentzVector *candpKpi=0x0;
 524    TLorentzVector *candpiKp=0x0;
 525
 526    invmasspKpi=InvMass(p1,p2,p3,candpKpi);
 527    invmasspiKp=InvMass(p3,p2,p1,candpiKp);
 528
 529    if(invmasspKpi < fThreesigmas[0] || invmasspKpi > fThreesigmas[1] ) hyppKpi=kFALSE;
 530    if(invmasspiKp < fThreesigmas[0] || invmasspiKp > fThreesigmas[1] ) hyppiKp=kFALSE;
 531
 532    if(!hyppKpi && !hyppiKp) {
 533       //Printf("Out of 3 sigmas Lc");
 534       delete candpiKp;
 535       delete candpKpi;
 536
 537       return;
 538    }
 539    // if(hyppKpi) Printf("INV MASS pKpi %f",invmasspKpi);
 540    // if(hyppiKp) Printf("piKp inv mass %f",invmasspiKp);
 541    hrej->Fill(6);
 542    //Double_t pcandpKpi=candpKpi->P();
 543    //Double_t pcandpiKp=candpiKp->P();
 544    Double_t ptcandpKpi=candpKpi->Pt();
 545    Double_t ptcandpiKp=candpiKp->Pt();
 546    //Printf("Filling hstograms");
 547    //Fill histograms
 548    if(etamax<0.8){
 549       if(hyppKpi){
 550          hPtEtaCandUpg->Fill(ptcandpKpi);
 551          hMassEtaCandUpg->Fill(candpKpi->M());
 552          hPtEtaCandMid->Fill(ptcandpKpi);
 553          hMassEtaCandMid->Fill(candpKpi->M());
 554       }
 555       if(hyppiKp) {
 556          hPtEtaCandMid->Fill(ptcandpiKp);
 557          hMassEtaCandMid->Fill(candpiKp->M());
 558          hPtEtaCandUpg->Fill(ptcandpiKp);
 559          hMassEtaCandUpg->Fill(candpiKp->M());
 560       }
 561    } else{ //eta >0.8 && <fAbsEta
 562       if(hyppKpi){
 563          hPtEtaCandUpg->Fill(ptcandpKpi);
 564          hMassEtaCandUpg->Fill(candpKpi->M());
 565       }
 566       if(hyppiKp) {
 567          hPtEtaCandUpg->Fill(ptcandpiKp);
 568          hMassEtaCandUpg->Fill(candpiKp->M());
 569
 570       }
 571
 572    }
 573
 574    delete candpiKp;
 575    delete candpKpi;
 576
 577 }
 578
 579 //__________________________________________________________________
 580 void AliAnalysisTaskCountLcEta::Terminate(Option_t *) {
 581
 582
 583 }
 584
 585 void AliAnalysisTaskCountLcEta::FillHistosL(TParticle *part, AliStack* stack){
 586
 587   //Histograms
 588   TH1F* hPtLc3Prongs=(TH1F*)fOutList->FindObject("hLc3Prongs");
 589   TH1F* hPLc3Prongs=(TH1F*)fOutList->FindObject("hPLc3Prongs");
 590
 591   TH1F* hPtpi=(TH1F*)fOutList->FindObject("hPtpi");
 592   TH1F* hPpi=(TH1F*)fOutList->FindObject("hPpi");
 593   TH1F* hPtK=(TH1F*)fOutList->FindObject("hPtK");
 594   TH1F* hPK=(TH1F*)fOutList->FindObject("hPK");
 595   TH1F* hPtp=(TH1F*)fOutList->FindObject("hPtp");
 596   TH1F* hPp=(TH1F*)fOutList->FindObject("hPp");
 597   TH1F* hLc3ProngsInEta=(TH1F*)fOutList->FindObject("hLc3ProngsInEta");
 598   TH1F* hLc3DaughInEta=(TH1F*)fOutList->FindObject("hLc3DaughInEta");
 599   TH1F* hLc2DaughInEta=(TH1F*)fOutList->FindObject("hLc2DaughInEta");
 600   TH1F* hLc1DaughInEta=(TH1F*)fOutList->FindObject("hLc1DaughInEta");
 601
 602   TH1F* hEtaLc=(TH1F*)fOutList->FindObject("hEtaLc");
 603   TH1F* hYLc=(TH1F*)fOutList->FindObject("hYLc");
 604
 605   TH1F* hLcpiInEta=(TH1F*)fOutList->FindObject("hLcpiInEta");
 606   TH1F* hLcpInEta=(TH1F*)fOutList->FindObject("hLcpInEta");
 607   TH1F* hLcKInEta=(TH1F*)fOutList->FindObject("hLcKInEta");
 608   TH1F* hLcpiKInEta=(TH1F*)fOutList->FindObject("hLcpiKInEta");
 609   TH1F* hLcpipInEta=(TH1F*)fOutList->FindObject("hLcpipInEta");
 610   TH1F* hLcKpInEta=(TH1F*)fOutList->FindObject("hLcKpInEta");
 611   TH1F* hLcpKpiInEta=(TH1F*)fOutList->FindObject("hLcpKpiInEta");
 612
 613
 614
 615   Double_t pt_part=part->Pt();
 616   Double_t p_part=part->P();
 617   Double_t eta_part=part->Eta();
 618   Double_t y_part=part->Y();
 619   Int_t nDaugh = part->GetNDaughters();
 620   //printf("Fillhistos L\n");
 621   if(nDaugh<2) return;
 622   if(nDaugh>3) return;
 623   //Printf("Lc in 3 prongs");
 624   TParticle* pdaugh1 = stack->Particle(part->GetFirstDaughter());
 625   if(!pdaugh1) return;
 626   Int_t number1 = TMath::Abs(pdaugh1->GetPdgCode());
 627   TParticle* pdaugh2 = stack->Particle(part->GetLastDaughter());
 628   if(!pdaugh2) return;
 629   Int_t number2 = TMath::Abs(pdaugh2->GetPdgCode());
 630
 631   Double_t eta1=0.;
 632   Double_t eta2=0.;
 633   Double_t eta3=0.;
 634   Bool_t okLambdac=kFALSE;
 635   Int_t pdgs[3]={0,0,0};
 636
 637   //TMath::Abs(pdaugh2->Eta());
 638
 639   //Lambda in 3 prongs
 640
 641   Double_t mom[3];
 642   Double_t mom_t[3];
 643   if(nDaugh==3){
 644     //Printf("Pt part %f",pt_part);
 645     Int_t thirdDaugh=part->GetLastDaughter()-1;
 646     TParticle* pdaugh3 = stack->Particle(thirdDaugh);
 647     //printf("Fillhistos L 3 daugh\n");
 648     if(!pdaugh3) return;
 649     mom[0]=pdaugh1->P();
 650     mom[1]=pdaugh2->P();
 651     mom[2]=pdaugh3->P();
 652     mom_t[0]=pdaugh1->Pt();
 653     mom_t[1]=pdaugh2->Pt();
 654     mom_t[2]=pdaugh3->Pt();
 655     eta1=TMath::Abs(pdaugh1->Eta());
 656     eta2=TMath::Abs(pdaugh2->Eta());
 657     eta3=TMath::Abs(pdaugh3->Eta());
 658     Int_t number3 = TMath::Abs(pdaugh3->GetPdgCode());
 659     pdgs[0]=number1; pdgs[1]=number2; pdgs[2]=number3;
 660     okLambdac=kFALSE;
 661     if(number1==211 && number2==321 && number3==2212) okLambdac=kTRUE;
 662     if(number1==321 && number2==211 && number3==2212) okLambdac=kTRUE;
 663     if(number1==321 && number2==2212 && number3==211) okLambdac=kTRUE;
 664     if(number1==211 && number2==2212 && number3==321) okLambdac=kTRUE;
 665     if(number1==2212 && number2==211 && number3==321) okLambdac=kTRUE;
 666     if(number1==2212 && number2==321 && number3==211) okLambdac=kTRUE;
 667   }
 668   if(nDaugh==2){
 669
 670     //Lambda resonant
 671     //Lambda -> p K*0
 672     Int_t nfiglieK=0;
 673     if(number1==2212 && number2==313){
 674       nfiglieK=pdaugh2->GetNDaughters();
 675       if(nfiglieK!=2) return;
 676       TParticle* pdaughK1 = stack->Particle(pdaugh2->GetFirstDaughter());
 677       TParticle* pdaughK2 = stack->Particle(pdaugh2->GetLastDaughter());
 678       if(!pdaughK1) return;
 679       if(!pdaughK2) return;
 680       Int_t number2K=TMath::Abs(pdaughK1->GetPdgCode());
 681       Int_t number3K=TMath::Abs(pdaughK2->GetPdgCode());
 682       mom[0]=pdaugh1->P();
 683       mom[1]=pdaughK1->P();
 684       mom[2]=pdaughK2->P();
 685       mom_t[0]=pdaugh1->Pt();
 686       mom_t[1]=pdaughK1->Pt();
 687       mom_t[2]=pdaughK2->Pt();
 688       pdgs[0]=number1; pdgs[1]=number2K; pdgs[2]=number3K;
 689       okLambdac=kFALSE;
 690       if(number1==211 && number2K==321 && number3K==2212) okLambdac=kTRUE;
 691       if(number1==321 && number2K==211 && number3K==2212) okLambdac=kTRUE;
 692       if(number1==321 && number2K==2212 && number3K==211) okLambdac=kTRUE;
 693       if(number1==211 && number2K==2212 && number3K==321) okLambdac=kTRUE;
 694       if(number1==2212 && number2K==211 && number3K==321) okLambdac=kTRUE;
 695       if(number1==2212 && number2K==321 && number3K==211) okLambdac=kTRUE;
 696       eta1=TMath::Abs(pdaugh1->Eta());
 697       eta2=TMath::Abs(pdaughK1->Eta());
 698       eta3=TMath::Abs(pdaughK2->Eta());
 699     }
 700
 701     if(number1==313 && number2==2212){
 702       nfiglieK=pdaugh1->GetNDaughters();
 703       if(nfiglieK!=2) return;
 704       TParticle* pdaughK1 = stack->Particle(pdaugh1->GetFirstDaughter());
 705       TParticle* pdaughK2 = stack->Particle(pdaugh1->GetLastDaughter());
 706       if(!pdaughK1) return;
 707       if(!pdaughK2) return;
 708       Int_t number2K=TMath::Abs(pdaughK1->GetPdgCode());
 709       Int_t number3K=TMath::Abs(pdaughK2->GetPdgCode());
 710       mom[0]=pdaugh2->P();
 711       mom[1]=pdaughK1->P();
 712       mom[2]=pdaughK2->P();
 713       mom_t[0]=pdaugh2->Pt();
 714       mom_t[1]=pdaughK1->Pt();
 715       mom_t[2]=pdaughK2->Pt();
 716       pdgs[0]=number2; pdgs[1]=number2K; pdgs[2]=number3K;
 717       okLambdac=kFALSE;
 718       if(number2==211 && number2K==321 && number3K==2212) okLambdac=kTRUE;
 719       if(number2==321 && number2K==211 && number3K==2212) okLambdac=kTRUE;
 720       if(number2==321 && number2K==2212 && number3K==211) okLambdac=kTRUE;
 721       if(number2==211 && number2K==2212 && number3K==321) okLambdac=kTRUE;
 722       if(number2==2212 && number2K==211 && number3K==321) okLambdac=kTRUE;
 723       if(number2==2212 && number2K==321 && number3K==211) okLambdac=kTRUE;
 724       eta1=TMath::Abs(pdaugh2->Eta());
 725       eta2=TMath::Abs(pdaughK1->Eta());
 726       eta3=TMath::Abs(pdaughK2->Eta());
 727     }
 728     //Lambda -> Delta++ k
 729     Int_t nfiglieDelta=0;
 730     if(number1==321 && number2==2224){
 731       nfiglieDelta=pdaugh2->GetNDaughters();
 732       if(nfiglieDelta!=2) return;
 733       TParticle *pdaughD1=stack->Particle(pdaugh2->GetFirstDaughter());
 734       TParticle *pdaughD2=stack->Particle(pdaugh2->GetLastDaughter());
 735       if(!pdaughD1) return;
 736       if(!pdaughD2) return;
 737       Int_t number2D=TMath::Abs(pdaughD1->GetPdgCode());
 738       Int_t number3D=TMath::Abs(pdaughD2->GetPdgCode());
 739       mom[0]=pdaugh1->P();
 740       mom[1]=pdaughD1->P();
 741       mom[2]=pdaughD2->P();
 742       mom_t[0]=pdaugh1->Pt();
 743       mom_t[1]=pdaughD1->Pt();
 744       mom_t[2]=pdaughD2->Pt();
 745       okLambdac=kFALSE;
 746       pdgs[0]=number1; pdgs[1]=number2D; pdgs[2]=number3D;
 747       if(number1==211 && number2D==321 && number3D==2212) okLambdac=kTRUE;
 748       if(number1==321 && number2D==211 && number3D==2212) okLambdac=kTRUE;
 749       if(number1==321 && number2D==2212 && number3D==211) okLambdac=kTRUE;
 750       if(number1==211 && number2D==2212 && number3D==321) okLambdac=kTRUE;
 751       if(number1==2212 && number2D==211 && number3D==321) okLambdac=kTRUE;
 752       if(number1==2212 && number2D==321 && number3D==211) okLambdac=kTRUE;
 753       eta1=TMath::Abs(pdaugh1->Eta());
 754       eta2=TMath::Abs(pdaughD1->Eta());
 755       eta3=TMath::Abs(pdaughD2->Eta());
 756     }
 757     if(number1==2224 && number2==321){
 758       nfiglieDelta=pdaugh1->GetNDaughters();
 759       if(nfiglieDelta!=2) return;
 760       TParticle* pdaughD1 = stack->Particle(pdaugh1->GetFirstDaughter());
 761       TParticle* pdaughD2 = stack->Particle(pdaugh1->GetLastDaughter());
 762       if(!pdaughD1) return;
 763       if(!pdaughD2) return;
 764       Int_t number2D=TMath::Abs(pdaughD1->GetPdgCode());
 765       Int_t number3D=TMath::Abs(pdaughD2->GetPdgCode());
 766       mom[0]=pdaugh2->P();
 767       mom[1]=pdaughD1->P();
 768       mom[2]=pdaughD2->P();
 769       mom_t[0]=pdaugh2->Pt();
 770       mom_t[1]=pdaughD1->Pt();
 771       mom_t[2]=pdaughD2->Pt();
 772       pdgs[0]=number2; pdgs[1]=number2D; pdgs[2]=number3D;
 773       okLambdac=kFALSE;
 774       if(number2==211 && number2D==321 && number3D==2212) okLambdac=kTRUE;
 775       if(number2==321 && number2D==211 && number3D==2212) okLambdac=kTRUE;
 776       if(number2==321 && number2D==2212 && number3D==211) okLambdac=kTRUE;
 777       if(number2==211 && number2D==2212 && number3D==321) okLambdac=kTRUE;
 778       if(number2==2212 && number2D==211 && number3D==321) okLambdac=kTRUE;
 779       if(number2==2212 && number2D==321 && number3D==211) okLambdac=kTRUE;
 780       eta1=TMath::Abs(pdaugh2->Eta());
 781       eta2=TMath::Abs(pdaughD1->Eta());
 782       eta3=TMath::Abs(pdaughD2->Eta());
 783     }
 784
 785     //Lambdac -> Lambda(1520) pi
 786     Int_t nfiglieLa=0;
 787     if(number1==3124 && number2==211){
 788       nfiglieLa=pdaugh1->GetNDaughters();
 789       if(nfiglieLa!=2) return;
 790       TParticle *pdaughL1=stack->Particle(pdaugh1->GetFirstDaughter());
 791       TParticle *pdaughL2=stack->Particle(pdaugh1->GetLastDaughter());
 792       if(!pdaughL1) return;
 793       if(!pdaughL2) return;
 794       Int_t number2L=TMath::Abs(pdaughL1->GetPdgCode());
 795       Int_t number3L=TMath::Abs(pdaughL2->GetPdgCode());
 796       mom[0]=pdaugh2->P();
 797       mom[1]=pdaughL1->P();
 798       mom[2]=pdaughL2->P();
 799       mom_t[0]=pdaugh2->Pt();
 800       mom_t[1]=pdaughL1->Pt();
 801       mom_t[2]=pdaughL2->Pt();
 802       pdgs[0]=number2; pdgs[1]=number2L; pdgs[2]=number3L;
 803       okLambdac=kFALSE;
 804       if(number2==211 && number2L==321 && number3L==2212) okLambdac=kTRUE;
 805       if(number2==321 && number2L==211 && number3L==2212) okLambdac=kTRUE;
 806       if(number2==321 && number2L==2212 && number3L==211) okLambdac=kTRUE;
 807       if(number2==211 && number2L==2212 && number3L==321) okLambdac=kTRUE;
 808       if(number2==2212 && number2L==211 && number3L==321) okLambdac=kTRUE;
 809       if(number2==2212 && number2L==321 && number3L==211) okLambdac=kTRUE;
 810       eta1=TMath::Abs(pdaugh2->Eta());
 811       eta2=TMath::Abs(pdaughL1->Eta());
 812       eta3=TMath::Abs(pdaughL2->Eta());
 813     }
 814     if(number1==211 && number2==3124){
 815       nfiglieLa=pdaugh2->GetNDaughters();
 816       if(nfiglieLa!=2) return;
 817       TParticle *pdaughL1=stack->Particle(pdaugh2->GetFirstDaughter());
 818       TParticle *pdaughL2=stack->Particle(pdaugh2->GetLastDaughter());
 819       if(!pdaughL1) return;
 820       if(!pdaughL2) return;
 821       Int_t number2L=TMath::Abs(pdaughL1->GetPdgCode());
 822       Int_t number3L=TMath::Abs(pdaughL2->GetPdgCode());
 823       mom[0]=pdaugh1->P();
 824       mom[1]=pdaughL1->P();
 825       mom[2]=pdaughL2->P();
 826       mom_t[0]=pdaugh1->Pt();
 827       mom_t[1]=pdaughL1->Pt();
 828       mom_t[2]=pdaughL2->Pt();
 829       pdgs[0]=number1; pdgs[1]=number2L; pdgs[2]=number3L;
 830       okLambdac=kFALSE;
 831       if(number1==211 && number2L==321 && number3L==2212) okLambdac=kTRUE;
 832       if(number1==321 && number2L==211 && number3L==2212) okLambdac=kTRUE;
 833       if(number1==321 && number2L==2212 && number3L==211) okLambdac=kTRUE;
 834       if(number1==211 && number2L==2212 && number3L==321) okLambdac=kTRUE;
 835       if(number1==2212 && number2L==211 && number3L==321) okLambdac=kTRUE;
 836       if(number1==2212 && number2L==321 && number3L==211) okLambdac=kTRUE;
 837       eta1=TMath::Abs(pdaugh1->Eta());
 838       eta2=TMath::Abs(pdaughL1->Eta());
 839       eta3=TMath::Abs(pdaughL2->Eta());
 840     }
 841
 842   }
 843   if(!okLambdac) return;
 844   //Printf("Filling");
 845   hPtLc3Prongs->Fill(pt_part);
 846   hPLc3Prongs->Fill(p_part);
 847   hEtaLc->Fill(eta_part);
 848   hYLc->Fill(y_part);
 849   if(eta_part<fEtaAbs) hLc3ProngsInEta->Fill(pt_part);
 850   for(Int_t iind=0;iind<3;iind++){
 851     if(pdgs[iind]==211) {
 852       hPtpi->Fill(mom_t[iind]);
 853       hPpi->Fill(mom[iind]);
 854     }
 855     if(pdgs[iind]==321) {
 856       hPtK->Fill(mom_t[iind]);
 857       hPK->Fill(mom[iind]);
 858     }
 859     if(pdgs[iind]==2212) {
 860       hPtp->Fill(mom_t[iind]);
 861       hPp->Fill(mom[iind]);
 862     }
 863   }
 864
 865   if(eta1<fEtaAbs && eta2<fEtaAbs && eta3<fEtaAbs){
 866     hLc3DaughInEta->Fill(pt_part);
 867   }
 868
 869   if((eta1<fEtaAbs && eta2<fEtaAbs &&  eta3>fEtaAbs) || (eta1<fEtaAbs && eta3<fEtaAbs &&  eta2>fEtaAbs) || (eta3<fEtaAbs && eta2<fEtaAbs &&  eta1>fEtaAbs)){
 870     hLc2DaughInEta->Fill(pt_part);
 871
 872     if(eta1<fEtaAbsMax && eta2<fEtaAbsMax && eta3<fEtaAbsMax){ //upper limit in eta (default 1.5)
 873       //PID : 1 id particle out of acceptance
 874       Int_t id=-1;
 875       if(eta1>fEtaAbs) id=0;
 876       if(eta2>fEtaAbs) id=1;
 877       if(eta3>fEtaAbs) id=2;
 878       if(pdgs[id]==211) hLcpiInEta->Fill(pt_part);
 879       if(pdgs[id]==321) hLcKInEta->Fill(pt_part);
 880       if(pdgs[id]==2212) hLcpInEta->Fill(pt_part);
 881     }
 882   }
 883
 884   if( (eta1<fEtaAbs && eta2>fEtaAbs && eta3>fEtaAbs) || (eta2<fEtaAbs && eta1>fEtaAbs && eta3>fEtaAbs) || (eta3<fEtaAbs && eta2>fEtaAbs &&  eta1>fEtaAbs)){
 885     hLc1DaughInEta->Fill(pt_part);
 886
 887     if(eta1<fEtaAbsMax && eta2<fEtaAbsMax && eta3<fEtaAbsMax){ //upper limit in eta (default 1.5)
 888       //PID : 2 id particles out of acceptance
 889       Int_t id1=-1,id2=-1;
 890       if(eta1>fEtaAbs && eta2>fEtaAbs) {
 891         id1=0;
 892         id2=1;
 893       }
 894       if(eta2>fEtaAbs && eta3>fEtaAbs) {
 895         id1=1;
 896         id2=2;
 897       }
 898       if(eta1>fEtaAbs && eta3>fEtaAbs) {
 899         id1=0;
 900         id2=2;
 901       }
 902       if((pdgs[id1]==211 && pdgs[id2]==321) || (pdgs[id1]==321 && pdgs[id2]==211))
 903         hLcpiKInEta->Fill(pt_part);
 904       if((pdgs[id1]==321 && pdgs[id2]==2212) || (pdgs[id2]==321 && pdgs[id1]==2212))
 905         hLcKpInEta->Fill(pt_part);
 906       if((pdgs[id1]==211 && pdgs[id2]==2212) || (pdgs[id2]==211 && pdgs[id1]==2212) )
 907         hLcpipInEta->Fill(pt_part);
 908     }
 909   }
 910
 911   if(eta1<fEtaAbsMax && eta2<fEtaAbsMax && eta3<fEtaAbsMax){ //upper limit in eta (default 1.5)
 912     //PID : 3 id particles out of acceptance
 913     if (eta1>fEtaAbs && eta2>fEtaAbs && eta3>fEtaAbs) {
 914       hLcpKpiInEta->Fill(pt_part);
 915     }
 916   }
 917
 918
 919
 920   return;
 921 }
 922
 923