AliPhysics  a9863a5 (a9863a5)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskSEDplus.cxx
Go to the documentation of this file.
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 // Class AliAnalysisTaskSEDplus
20 // AliAnalysisTaskSE for the D+ candidates Invariant Mass Histogram and
21 // comparison of heavy-flavour decay candidates
22 // to MC truth (kinematics stored in the AOD)
23 // Authors: Renu Bala, bala@to.infn.it
24 // F. Prino, prino@to.infn.it
25 // G. Ortona, ortona@to.infn.it
27 
28 #include <TClonesArray.h>
29 #include <TNtuple.h>
30 #include <TCanvas.h>
31 #include <TList.h>
32 #include <TString.h>
33 #include <TDatabasePDG.h>
34 
35 #include "AliAnalysisManager.h"
37 #include "AliAODHandler.h"
38 #include "AliAODEvent.h"
39 #include "AliAODVertex.h"
40 #include "AliAODTrack.h"
42 #include "AliAnalysisVertexingHF.h"
43 #include "AliAnalysisTaskSE.h"
44 #include "AliAnalysisTaskSEDplus.h"
46 #include "AliVertexingHFUtils.h"
47 
51 
52 //________________________________________________________________________
54 AliAnalysisTaskSE(),
55  fOutput(0),
56  fHistNEvents(0),
57  fMCAccPrompt(0),
58  fMCAccBFeed(0),
59  fPtVsMassNoPid(0),
60  fPtVsMass(0),
61  fYVsPtNoPid(0),
62  fYVsPt(0),
63  fYVsPtSigNoPid(0),
64  fYVsPtSig(0),
65  fPhiEtaCand(0),
66  fPhiEtaCandSigReg(0),
67  fSPDMult(0),
68  fNtupleDplus(0),
69  fUpmasslimit(1.965),
70  fLowmasslimit(1.765),
71  fNPtBins(0),
72  fBinWidth(0.002),
73  fListCuts(0),
74  fRDCutsAnalysis(0),
75  fCounter(0),
76  fFillNtuple(0),
77  fReadMC(kFALSE),
78  fUseStrangeness(kFALSE),
79  fUseBit(kTRUE),
80  fCutsDistr(kFALSE),
81  fDoImpPar(kFALSE),
82  fStepMCAcc(kFALSE),
83  fUseQuarkTagInKine(kTRUE),
84  fNImpParBins(400),
85  fLowerImpPar(-1000.),
86  fHigherImpPar(1000.),
87  fDoLS(0),
88  fEtaSelection(0),
89  fSystem(0)
90 {
92 
93  for(Int_t i=0;i<3;i++){
94  fHistCentrality[i]=0;
95  fCorreld0Kd0pi[i]=0;
96  }
97 
98  for(Int_t i=0; i<5; i++)fHistMassPtImpPar[i]=0;
99 
100 
101  for(Int_t i=0;i<3*kMaxPtBins;i++){
102  fMassHistNoPid[i]=0;
103  fCosPHist[i]=0;
104  fDLenHist[i]=0;
105  fSumd02Hist[i]=0;
106  fSigVertHist[i]=0;
107  fPtMaxHist[i]=0;
108  fPtKHist[i]=0;
109  fPtpi1Hist[i]=0;
110  fPtpi2Hist[i]=0;
111  fDCAHist[i]=0;
112  fMassHist[i]=0;
113  fMassHistPlus[i]=0;
114  fMassHistMinus[i]=0;
115 
116  fDLxy[i]=0;
117  fCosxy[i]=0;
118  fCosPHistLS[i]=0;
119  fDLenHistLS[i]=0;
120  fSumd02HistLS[i]=0;
121  fSigVertHistLS[i]=0;
122  fPtMaxHistLS[i]=0;
123  fDCAHistLS[i]=0;
124  }
125  for(Int_t i=0;i<5*kMaxPtBins;i++){
126  fMassHistLS[i]=0;
127  }
128  for(Int_t i=0;i<kMaxPtBins+1;i++){
129  fArrayBinLimits[i]=0;
130  }
131 
132 }
133 
134 //________________________________________________________________________
135 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,AliRDHFCutsDplustoKpipi *dpluscutsana,Int_t fillNtuple):
136  AliAnalysisTaskSE(name),
137  fOutput(0),
138  fHistNEvents(0),
139  fMCAccPrompt(0),
140  fMCAccBFeed(0),
141  fPtVsMassNoPid(0),
142  fPtVsMass(0),
143  fYVsPtNoPid(0),
144  fYVsPt(0),
145  fYVsPtSigNoPid(0),
146  fYVsPtSig(0),
147  fPhiEtaCand(0),
148  fPhiEtaCandSigReg(0),
149  fSPDMult(0),
150  fNtupleDplus(0),
151  fUpmasslimit(1.965),
152  fLowmasslimit(1.765),
153  fNPtBins(0),
154  fBinWidth(0.002),
155  fListCuts(0),
156  fRDCutsAnalysis(dpluscutsana),
157  fCounter(0),
158  fFillNtuple(fillNtuple),
159  fReadMC(kFALSE),
160  fUseStrangeness(kFALSE),
161  fUseBit(kTRUE),
162  fCutsDistr(kFALSE),
163  fDoImpPar(kFALSE),
164  fStepMCAcc(kFALSE),
165  fUseQuarkTagInKine(kTRUE),
166  fNImpParBins(400),
167  fLowerImpPar(-1000.),
168  fHigherImpPar(1000.),
169  fDoLS(0),
170  fEtaSelection(0),
171  fSystem(0)
172 {
173  //
175  //
177 
178  for(Int_t i=0;i<3;i++){
179  fHistCentrality[i]=0;
180  fCorreld0Kd0pi[i]=0;
181  }
182 
183  for(Int_t i=0; i<5; i++)fHistMassPtImpPar[i]=0;
184 
185  for(Int_t i=0;i<3*kMaxPtBins;i++){
186  fMassHistNoPid[i]=0;
187  fCosPHist[i]=0;
188  fDLenHist[i]=0;
189  fSumd02Hist[i]=0;
190  fSigVertHist[i]=0;
191  fPtMaxHist[i]=0;
192  fPtKHist[i]=0;
193  fPtpi1Hist[i]=0;
194  fPtpi2Hist[i]=0;
195  fDCAHist[i]=0;
196  fMassHist[i]=0;
197  fMassHistPlus[i]=0;
198  fMassHistMinus[i]=0;
199 
200  fDLxy[i]=0;
201  fCosxy[i]=0;
202  fCosPHistLS[i]=0;
203  fDLenHistLS[i]=0;
204  fSumd02HistLS[i]=0;
205  fSigVertHistLS[i]=0;
206  fPtMaxHistLS[i]=0;
207  fDCAHistLS[i]=0;
208  }
209  for(Int_t i=0;i<5*kMaxPtBins;i++){
210  fMassHistLS[i]=0;
211  }
212  for(Int_t i=0;i<kMaxPtBins+1;i++){
213  fArrayBinLimits[i]=0;
214  }
215 
216 
217  // Default constructor
218  // Output slot #1 writes into a TList container
219  DefineOutput(1,TList::Class()); //My private output
220  // Output slot #2 writes cut to private output
221  // DefineOutput(2,AliRDHFCutsDplustoKpipi::Class());
222  DefineOutput(2,TList::Class());
223  // Output slot #3 writes cut to private output
224  DefineOutput(3,AliNormalizationCounter::Class());
225 
226  if(fFillNtuple){
227  // Output slot #4 writes into a TNtuple container
228  DefineOutput(4,TNtuple::Class()); //My private output
229  }
230 }
231 
232 //________________________________________________________________________
234 {
235  //
237  //
238  if(fOutput && !fOutput->IsOwner()){
239  delete fHistNEvents;
240  for(Int_t i=0;i<3*fNPtBins;i++){
241  delete fMassHistNoPid[i];
242  delete fCosPHist[i];
243  delete fDLenHist[i];
244  delete fSumd02Hist[i];
245  delete fSigVertHist[i];
246  delete fPtMaxHist[i];
247  delete fPtKHist[i];
248  delete fPtpi1Hist[i];
249  delete fPtpi2Hist[i];
250  delete fDCAHist[i];
251  delete fDLxy[i];
252  delete fCosxy[i];
253  delete fMassHist[i];
254  delete fMassHistPlus[i];
255  delete fMassHistMinus[i];
256  delete fCosPHistLS[i];
257  delete fDLenHistLS[i];
258  delete fSumd02HistLS[i];
259  delete fSigVertHistLS[i];
260  delete fPtMaxHistLS[i];
261  delete fDCAHistLS[i];
262  }
263  for(Int_t i=0;i<5*fNPtBins;i++){
264  delete fMassHistLS[i];
265  }
266  for(Int_t i=0;i<3;i++){
267  delete fCorreld0Kd0pi[i];
268  delete fHistCentrality[i];
269  }
270  for(Int_t i=0;i<5;i++){
271  delete fHistMassPtImpPar[i];
272  }
273  delete fPtVsMassNoPid;
274  delete fPtVsMass;
275  delete fYVsPtNoPid;
276  delete fYVsPt;
277  delete fYVsPtSigNoPid;
278  delete fYVsPtSig;
279  delete fPhiEtaCand;
280  delete fPhiEtaCandSigReg;
281  delete fSPDMult;
282  delete fMCAccPrompt;
283  delete fMCAccBFeed;
284  }
285 
286  delete fOutput;
287  delete fNtupleDplus;
288  delete fListCuts;
289  delete fRDCutsAnalysis;
290  delete fCounter;
291 
292 }
293 //_________________________________________________________________
296  Float_t bw=GetBinWidth();
297  fUpmasslimit = 1.865+range;
298  fLowmasslimit = 1.865-range;
299  SetBinWidth(bw);
300 }
301 //_________________________________________________________________
302 void AliAnalysisTaskSEDplus::SetMassLimits(Float_t lowlimit, Float_t uplimit){
304  if(uplimit>lowlimit)
305  {
306  Float_t bw=GetBinWidth();
307  fUpmasslimit = lowlimit;
308  fLowmasslimit = uplimit;
309  SetBinWidth(bw);
310  }
311 }
312 //________________________________________________________________
314  Float_t width=w;
315  Int_t nbins=(Int_t)((fUpmasslimit-fLowmasslimit)/width+0.5);
316  Int_t missingbins=4-nbins%4;
317  nbins=nbins+missingbins;
318  width=(fUpmasslimit-fLowmasslimit)/nbins;
319  if(missingbins!=0){
320  printf("AliAnalysisTaskSEDplus::SetBinWidth: W-bin width of %f will produce histograms not rebinnable by 4. New width set to %f\n",w,width);
321  }
322  else{
323  if(fDebug>1) printf("AliAnalysisTaskSEDplus::SetBinWidth: width set to %f\n",width);
324  }
325  fBinWidth=width;
326 }
327 //_________________________________________________________________
329  return (Int_t)((fUpmasslimit-fLowmasslimit)/fBinWidth+0.5);
330 }
331 //_________________________________________________________________
332 void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS){
333  //
334  //
336  //
337  if(fDebug>1)printf("started LS\n");
338 
339  //histograms for like sign
340  Int_t nbins=GetNBinsHistos();;
341  TH1F *histLSPlus = new TH1F("LSPlus","LSPlus",nbins,fLowmasslimit,fUpmasslimit);
342  TH1F *histLSMinus = new TH1F("LSMinus","LSMinus",nbins,fLowmasslimit,fUpmasslimit);
343 
344  Int_t nPosTrks=0,nNegTrks=0;
345  Int_t nOStriplets = arrayOppositeSign->GetEntriesFast();
346  Int_t nDplusLS=0;
347  Int_t nDminusLS=0;
348  Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
349  Int_t index=0;
350 
351  // loop over like sign candidates
352  for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
353  AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
354  if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts))continue;
355  Bool_t unsetvtx=kFALSE;
356  if(!d->GetOwnPrimaryVtx()) {
357  d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
358  unsetvtx=kTRUE;
359  }
360  Double_t ptCand = d->Pt();
361  Int_t iPtBin = fRDCutsAnalysis->PtBin(ptCand);
362  if(iPtBin<0)continue;
363  Int_t sign= d->GetCharge();
364  if(sign>0){
365  nPosTrks++;
366  }else{
367  nNegTrks++;
368  }
370  Int_t passTightCuts=fRDCutsAnalysis->GetIsSelectedCuts();
371  if(passTightCuts>0){
372  Float_t invMass = d->InvMassDplus();
373  index=GetLSHistoIndex(iPtBin);
374  fMassHistLS[index+1]->Fill(invMass);
375  if(sign>0){
376  histLSPlus->Fill(invMass);
377  nDplusLS++;
378  }else{
379  histLSMinus->Fill(invMass);
380  nDminusLS++;
381  }
382  if(fCutsDistr){
383  Double_t dlen=d->DecayLength();
384  Double_t cosp=d->CosPointingAngle();
385  Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
386  Double_t dca=d->GetDCA();
387  Double_t sigvert=d->GetSigmaVert();
388  Double_t ptmax=0;
389  for(Int_t i=0;i<3;i++){
390  if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
391  }
392  fCosPHistLS[iPtBin]->Fill(cosp);
393  fDLenHistLS[iPtBin]->Fill(dlen);
394  fSumd02HistLS[iPtBin]->Fill(sumD02);
395  fSigVertHistLS[iPtBin]->Fill(sigvert);
396  fPtMaxHistLS[iPtBin]->Fill(ptmax);
397  fDCAHistLS[iPtBin]->Fill(dca);
398  }
399  }
400  if(unsetvtx) d->UnsetOwnPrimaryVtx();
401  }
402  //wei:normalized to the number of combinations (production)
403  //wei2:normalized to the number of LS/OS (production)
404  //wei3:normalized to the number of LS/OS (analysis)
405  //wei4:normalized to the number of combinations (analysis)
406  Float_t wei2=0;
407  if(nLikeSign!=0)wei2 = (Float_t)nOStriplets/(Float_t)nLikeSign;
408  Float_t wei3=0;
409  if(nDplusLS!=0)wei3 = (Float_t)nDplusOS/(Float_t)(nDplusLS+nDminusLS);
410  Float_t weiplus=1.,weiminus=1.;
411  Float_t wei4plus=1.,wei4minus=1.;
412  //wei* should be automatically protected, since to get a triplet there must be at least 3 good tracks in the event
413  if(nPosTrks>2)weiplus=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.);
414  if(nDplusLS>2)wei4plus=3.*(Float_t)nDminusLS/((Float_t)nDplusLS-2.);
415  if(nNegTrks>2)weiminus=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.);
416  if(nDminusLS>2)wei4minus=3.*(Float_t)nDplusLS/((Float_t)nDminusLS-2.);
417 
418  fMassHistLS[index]->Add(histLSPlus,weiplus);
419  fMassHistLS[index]->Add(histLSMinus,weiminus);
420  fMassHistLS[index+2]->Add(histLSPlus,wei2);
421  fMassHistLS[index+2]->Add(histLSMinus,wei2);
422  fMassHistLS[index+3]->Add(histLSPlus,wei3);
423  fMassHistLS[index+3]->Add(histLSMinus,wei3);
424  fMassHistLS[index+4]->Add(histLSPlus,wei4plus);
425  fMassHistLS[index+4]->Add(histLSMinus,wei4minus);
426 
427  delete histLSPlus;histLSPlus=0;
428  delete histLSMinus;histLSMinus=0;
429 
430  if(fDebug>1) printf("LS analysis terminated\n");
431 }
432 
433 //__________________________________________
435  //
437  //
438  if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
439 
440  //PostData(2,fRDCutsloose);//we should then put those cuts in a tlist if we have more than 1
441  fListCuts=new TList();
442  fListCuts->SetOwner();
444  analysis->SetName("AnalysisCuts");
445 
446  fListCuts->Add(analysis);
447  PostData(2,fListCuts);
448 
449  return;
450 }
451 
452 //________________________________________________________________________
454 {
456  //
457  if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");
458 
459  // Several histograms are more conveniently managed in a TList
460  fOutput = new TList();
461  fOutput->SetOwner();
462  fOutput->SetName("OutputHistos");
463 
464  fHistNEvents = new TH1F("fHistNEvents", "number of events ",13,-0.5,12.5);
465  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
466  fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents accepted");
467  fHistNEvents->GetXaxis()->SetBinLabel(3,"Rejected due to trigger");
468  fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected pileup events");
469  fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to centrality");
470  fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vtxz");
471  fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to Physics Sel");
472  fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
473  fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
474  fHistNEvents->GetXaxis()->SetBinLabel(10,"D+ after topological cuts");
475  fHistNEvents->GetXaxis()->SetBinLabel(11,"D+ after Topological+PID cuts");
476  fHistNEvents->GetXaxis()->SetBinLabel(12,"D+ not on-the-fly reco");
477 
478  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
479  fHistNEvents->Sumw2();
480  fHistNEvents->SetMinimum(0);
481  fOutput->Add(fHistNEvents);
482 
483  TString hisname;
484  Int_t index=0;
485  Int_t nbins=GetNBinsHistos();
486  fHistCentrality[0]=new TH2F("hCentrMult","centrality",100,0.5,30000.5,40,0.,100.);
487  fHistCentrality[1]=new TH2F("hCentrMult(selectedCent)","centrality(selectedCent)",100,0.5,30000.5,40,0.,100.);
488  fHistCentrality[2]=new TH2F("hCentrMult(OutofCent)","centrality(OutofCent)",100,0.5,30000.5,40,0.,100.);
489  for(Int_t i=0;i<3;i++){
490  fHistCentrality[i]->Sumw2();
491  fOutput->Add(fHistCentrality[i]);
492  }
493  for(Int_t i=0;i<fNPtBins;i++){
494 
495  index=GetHistoIndex(i);
496 
497  hisname.Form("hMassNoPidPt%d",i);
498  fMassHistNoPid[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
499  fMassHistNoPid[index]->Sumw2();
500  hisname.Form("hCosPAllPt%d",i);
501  fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
502  fCosPHist[index]->Sumw2();
503  hisname.Form("hDLenAllPt%d",i);
504  fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
505  fDLenHist[index]->Sumw2();
506  hisname.Form("hSumd02AllPt%d",i);
507  fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
508  fSumd02Hist[index]->Sumw2();
509  hisname.Form("hSigVertAllPt%d",i);
510  fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
511  fSigVertHist[index]->Sumw2();
512  hisname.Form("hPtMaxAllPt%d",i);
513  fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
514  fPtMaxHist[index]->Sumw2();
515  hisname.Form("hPtKPt%d",i);
516  fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
517  fPtKHist[index]->Sumw2();
518  hisname.Form("hPtpi1Pt%d",i);
519  fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
520  fPtpi1Hist[index]->Sumw2();
521  hisname.Form("hPtpi2Pt%d",i);
522  fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
523  fPtpi2Hist[index]->Sumw2();
524  hisname.Form("hDCAAllPt%d",i);
525  fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
526  fDCAHist[index]->Sumw2();
527 
528  hisname.Form("hDLxyPt%d",i);
529  fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
530  fDLxy[index]->Sumw2();
531  hisname.Form("hCosxyPt%d",i);
532  fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
533  fCosxy[index]->Sumw2();
534 
535  hisname.Form("hMassPt%dTC",i);
536  fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
537  fMassHist[index]->Sumw2();
538  hisname.Form("hMassPt%dTCPlus",i);
539  fMassHistPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
540  fMassHistPlus[index]->Sumw2();
541  hisname.Form("hMassPt%dTCMinus",i);
542  fMassHistMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
543  fMassHistMinus[index]->Sumw2();
544 
545 
546 
547  index=GetSignalHistoIndex(i);
548  hisname.Form("hSigNoPidPt%d",i);
549  fMassHistNoPid[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
550  fMassHistNoPid[index]->Sumw2();
551  hisname.Form("hCosPSigPt%d",i);
552  fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
553  fCosPHist[index]->Sumw2();
554  hisname.Form("hDLenSigPt%d",i);
555  fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
556  fDLenHist[index]->Sumw2();
557  hisname.Form("hSumd02SigPt%d",i);
558  fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
559  fSumd02Hist[index]->Sumw2();
560  hisname.Form("hSigVertSigPt%d",i);
561  fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
562  fSigVertHist[index]->Sumw2();
563  hisname.Form("hPtMaxSigPt%d",i);
564  fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
565  fPtMaxHist[index]->Sumw2();
566  hisname.Form("hPtKSigPt%d",i);
567  fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
568  fPtKHist[index]->Sumw2();
569  hisname.Form("hPtpi1SigPt%d",i);
570  fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
571  fPtpi1Hist[index]->Sumw2();
572  hisname.Form("hPtpi2SigPt%d",i);
573  fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
574  fPtpi2Hist[index]->Sumw2();
575 
576  hisname.Form("hDCASigPt%d",i);
577  fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
578  fDCAHist[index]->Sumw2();
579 
580  hisname.Form("hDLxySigPt%d",i);
581  fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
582  fDLxy[index]->Sumw2();
583  hisname.Form("hCosxySigPt%d",i);
584  fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
585  fCosxy[index]->Sumw2();
586  hisname.Form("hSigPt%dTC",i);
587  fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
588  fMassHist[index]->Sumw2();
589  hisname.Form("hSigPt%dTCPlus",i);
590  fMassHistPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
591  fMassHistPlus[index]->Sumw2();
592  hisname.Form("hSigPt%dTCMinus",i);
593  fMassHistMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
594  fMassHistMinus[index]->Sumw2();
595 
596 
597  index=GetBackgroundHistoIndex(i);
598  hisname.Form("hBkgNoPidPt%d",i);
599  fMassHistNoPid[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
600  fMassHistNoPid[index]->Sumw2();
601  hisname.Form("hCosPBkgPt%d",i);
602  fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
603  fCosPHist[index]->Sumw2();
604  hisname.Form("hDLenBkgPt%d",i);
605  fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
606  fDLenHist[index]->Sumw2();
607  hisname.Form("hSumd02BkgPt%d",i);
608  fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
609  fSumd02Hist[index]->Sumw2();
610  hisname.Form("hSigVertBkgPt%d",i);
611  fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
612  fSigVertHist[index]->Sumw2();
613  hisname.Form("hPtMaxBkgPt%d",i);
614  fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
615  fPtMaxHist[index]->Sumw2();
616  hisname.Form("hPtKBkgPt%d",i);
617  fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
618  fPtKHist[index]->Sumw2();
619  hisname.Form("hPtpi1BkgPt%d",i);
620  fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
621  fPtpi1Hist[index]->Sumw2();
622  hisname.Form("hPtpi2BkgPt%d",i);
623  fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
624  fPtpi2Hist[index]->Sumw2();
625  hisname.Form("hDCABkgPt%d",i);
626  fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
627  fDCAHist[index]->Sumw2();
628 
629  hisname.Form("hDLxyBkgPt%d",i);
630  fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
631  fDLxy[index]->Sumw2();
632  hisname.Form("hCosxyBkgPt%d",i);
633  fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
634  fCosxy[index]->Sumw2();
635 
636 
637  hisname.Form("hBkgPt%dTC",i);
638  fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
639  fMassHist[index]->Sumw2();
640  hisname.Form("hBkgPt%dTCPlus",i);
641  fMassHistPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
642  fMassHistPlus[index]->Sumw2();
643  hisname.Form("hBkgPt%dTCMinus",i);
644  fMassHistMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
645  fMassHistMinus[index]->Sumw2();
646  }
647 
648 
649  for(Int_t i=0; i<3*fNPtBins; i++){
650  fOutput->Add(fMassHistNoPid[i]);
651  if(fCutsDistr){
652  fOutput->Add(fCosPHist[i]);
653  fOutput->Add(fDLenHist[i]);
654  fOutput->Add(fSumd02Hist[i]);
655  fOutput->Add(fSigVertHist[i]);
656  fOutput->Add(fPtMaxHist[i]);
657  fOutput->Add(fPtKHist[i]);
658  fOutput->Add(fPtpi1Hist[i]);
659  fOutput->Add(fPtpi2Hist[i]);
660  fOutput->Add(fDCAHist[i]);
661  fOutput->Add(fDLxy[i]);
662  fOutput->Add(fCosxy[i]);
663  }
664  fOutput->Add(fMassHist[i]);
665  fOutput->Add(fMassHistPlus[i]);
666  fOutput->Add(fMassHistMinus[i]);
667 
668  }
669 
670  fCorreld0Kd0pi[0]=new TH2F("hCorreld0Kd0piAll","",100,-0.02,0.02,100,-0.02,0.02);
671  fCorreld0Kd0pi[1]=new TH2F("hCorreld0Kd0piSig","",100,-0.02,0.02,100,-0.02,0.02);
672  fCorreld0Kd0pi[2]=new TH2F("hCorreld0Kd0piBkg","",100,-0.02,0.02,100,-0.02,0.02);
673  if(fCutsDistr){
674  for(Int_t i=0; i<3; i++){
675  fCorreld0Kd0pi[i]->Sumw2();
676  fOutput->Add(fCorreld0Kd0pi[i]);
677  }
678  }
679 
680 
681  fPtVsMassNoPid=new TH2F("hPtVsMassNoPid","PtVsMass (no PID)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
682  fPtVsMass=new TH2F("hPtVsMass","PtVsMass",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
683  fYVsPtNoPid=new TH3F("hYVsPtNoPid","YvsPt (no PID)",40,0.,20.,80,-2.,2.,nbins,fLowmasslimit,fUpmasslimit);
684  fYVsPt=new TH3F("hYVsPt","YvsPt",40,0.,20.,80,-2.,2.,nbins,fLowmasslimit,fUpmasslimit);
685  fYVsPtSigNoPid=new TH2F("hYVsPtSigNoPid","YvsPt (MC, only sig., no PID)",40,0.,20.,80,-2.,2.);
686  fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only Sig)",40,0.,20.,80,-2.,2.);
687  fPhiEtaCand=new TH2F("hPhiEtaCand","phi vs. eta candidates",20,-1.,1.,50,0.,2*TMath::Pi());
688  fPhiEtaCandSigReg=new TH2F("hPhiEtaCandSigReg","phi vs. eta candidates",20,-1.,1.,50,0.,2*TMath::Pi());
689 
690  Double_t maxmult;
691  if(fSystem==1) maxmult=5000;
692  else maxmult=200;
693  fSPDMult = new TH1F("hSPDMult", "Tracklets multiplicity; Tracklets ; Entries",200,0.,maxmult);
694  fOutput->Add(fPtVsMassNoPid);
695  fOutput->Add(fPtVsMass);
696  fOutput->Add(fYVsPtNoPid);
697  fOutput->Add(fYVsPt);
698  fOutput->Add(fYVsPtSigNoPid);
699  fOutput->Add(fYVsPtSig);
700  fOutput->Add(fPhiEtaCand);
702  fOutput->Add(fSPDMult);
703 
704 
705  //Counter for Normalization
706  TString normName="NormalizationCounter";
707  AliAnalysisDataContainer *cont = GetOutputSlot(3)->GetContainer();
708  if(cont)normName=(TString)cont->GetName();
709  fCounter = new AliNormalizationCounter(normName.Data());
710  fCounter->Init();
711 
715 
716  if(fFillNtuple==1){
717  OpenFile(4); // 4 is the slot number of the ntuple
718 
719 
720  fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:Pt:charge:piddau0:piddau1:piddau2:Ptpi:PtK:Ptpi2:mompi:momK:mompi2:cosp:cospxy:DecLen:NormDecLen:DecLenXY:NormDecLenXY:InvMass:sigvert:d0Pi:d0K:d0Pi2:maxdca:ntracks:centr:RunNumber:BadDau");
721 
722  }
723 
724  if(fFillNtuple==2){
725  OpenFile(4); // 4 is the slot number of the ntuple
726 
727 
728  fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Pt:InvMass:d0:origin");
729 
730  }
731 
732  return;
733 }
734 
735 //________________________________________________________________________
736 void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
737 {
740 
741  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
742 
743  TClonesArray *array3Prong = 0;
744  TClonesArray *arrayLikeSign =0;
745  if(!aod && AODEvent() && IsStandardAOD()) {
746  // In case there is an AOD handler writing a standard AOD, use the AOD
747  // event in memory rather than the input (ESD) event.
748  aod = dynamic_cast<AliAODEvent*> (AODEvent());
749  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
750  // have to taken from the AOD event hold by the AliAODExtension
751  AliAODHandler* aodHandler = (AliAODHandler*)
752  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
753  if(aodHandler->GetExtensions()) {
754  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
755  AliAODEvent *aodFromExt = ext->GetAOD();
756  array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
757  arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
758  }
759  } else if(aod) {
760  array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
761  arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
762  }
763 
764  if(!aod || !array3Prong) {
765  printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
766  return;
767  }
768  if(!arrayLikeSign && fDoLS) {
769  printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
770  return;
771  }
772 
773  // fix for temporary bug in ESDfilter
774  // the AODs with null vertex pointer didn't pass the PhysSel
775  if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
776 
777  //Store the event in AliNormalizationCounter->To disable for Pb-Pb? Add a flag to disable it?
779 
780  fHistNEvents->Fill(0); // count event
781  Int_t runNumber=aod->GetRunNumber();
782 
783  //Event selection
784  Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
785  Float_t ntracks=aod->GetNumberOfTracks();
786  Float_t evCentr=0;
788  fHistCentrality[0]->Fill(ntracks,evCentr);
789  if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(2);
790  if(fRDCutsAnalysis->GetWhyRejection()==1) fHistNEvents->Fill(3);
791  if(fRDCutsAnalysis->GetWhyRejection()==2){fHistNEvents->Fill(4);fHistCentrality[2]->Fill(ntracks,evCentr);}
792  if(fRDCutsAnalysis->GetWhyRejection()==6)fHistNEvents->Fill(5);
793  if(fRDCutsAnalysis->GetWhyRejection()==7)fHistNEvents->Fill(6);
794 
795  TClonesArray *arrayMC=0;
796  AliAODMCHeader *mcHeader=0;
797  // load MC particles
798  if(fReadMC){
799  arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
800  if(!arrayMC) {
801  printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
802  return;
803  }
804 
805  // load MC header
806  mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
807  if(!mcHeader) {
808  printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
809  return;
810  }
811  }
812  if(fReadMC && fStepMCAcc){
813  if(aod->GetTriggerMask()==0 &&
814  (runNumber>=195344 && runNumber<=195677)){
815  // protection for events with empty trigger mask in p-Pb
816  return;
817  }
819  // events not passing the centrality selection can be removed immediately.
820 
821  FillMCAcceptanceHistos(arrayMC, mcHeader);
822  }
823  // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
824  //TString trigclass=aod->GetFiredTriggerClasses();
825  // Post the data already here
826  PostData(1,fOutput);
827  if(!isEvSel)return;
828  Int_t tracklets=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
829  // printf("ntracklet===%d\n",tracklets);
830  fSPDMult->Fill(tracklets);
831 
832  fHistCentrality[1]->Fill(ntracks,evCentr);
833  fHistNEvents->Fill(1);
834 
835  // AOD primary vertex
836  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
837  // vtx1->Print();
838  // TString primTitle = vtx1->GetTitle();
839  //if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0)fHistNEvents->Fill(2);
840 
841 
842  Int_t n3Prong = array3Prong->GetEntriesFast();
843  // printf("Number of D+->Kpipi: %d and of tracks: %d\n",n3Prong,aod->GetNumberOfTracks());
844 
845  Int_t nOS=0;
846  Int_t index;
847  Int_t pdgDgDplustoKpipi[3]={321,211,211};
848 
849  // vHF object is needed to call the method that refills the missing info of the candidates
850  // if they have been deleted in dAOD reconstruction phase
851  // in order to reduce the size of the file
853  if(fDoLS>1){//Normalizations for LS
854  for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
855  AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
858  }
859  }
860  }else{//Standard analysis
861  // Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};//TO REMOVE
862  //Double_t *cutsDplus = new (Double_t*)fRDCuts->GetCuts();
863  Int_t nSelectednopid=0,nSelected=0;
864  for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
865  AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
866  fHistNEvents->Fill(7);
868  fHistNEvents->Fill(8);
869  continue;
870  }
871 
872  if(!(vHF->FillRecoCand(aod,d))) { //Fill the data members of the candidate only if they are empty.
873  fHistNEvents->Fill(12); //monitor how often this fails
874  continue;
875  }
876 
877  Int_t passTopolAndPIDCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
878 
879  if(!fRDCutsAnalysis->GetIsSelectedCuts()) continue;
880 
881  Double_t etaD=d->Eta();
882  Double_t phiD=d->Phi();
883  if(fEtaSelection!=0){
884  if(fEtaSelection==1 && etaD<0) continue;
885  if(fEtaSelection==-1 && etaD>0) continue;
886  }
887 
888  Bool_t unsetvtx=kFALSE;
889  if(!d->GetOwnPrimaryVtx()){
890  d->SetOwnPrimaryVtx(vtx1);
891  unsetvtx=kTRUE;
892  }
893 
894  Double_t ptCand = d->Pt();
895  Int_t iPtBin = fRDCutsAnalysis->PtBin(ptCand);
896 
897  Bool_t recVtx=kFALSE;
898  AliAODVertex *origownvtx=0x0;
900  if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
901  if(fRDCutsAnalysis->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
902  else fRDCutsAnalysis->CleanOwnPrimaryVtx(d,aod,origownvtx);
903  }
904 
905  Int_t labDp=-1;
906  Bool_t isPrimary=kFALSE;
907  Bool_t isFeeddown=kFALSE;
908  Float_t pdgCode=-2;
909  Float_t trueImpParXY=0.;
910  Double_t ptB=-1.5;
911  if(fReadMC){
912  labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
913  if(labDp>=0){
914  AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
915  Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,partDp,fUseQuarkTagInKine);//Prompt = 4, FeedDown = 5
916  pdgCode=TMath::Abs(partDp->GetPdgCode());
917  if(orig==4){
918  isPrimary=kTRUE;
919  isFeeddown=kFALSE;
920  }else if(orig==5){
921  isPrimary=kFALSE;
922  isFeeddown=kTRUE;
923  trueImpParXY=GetTrueImpactParameter(mcHeader,arrayMC,partDp)*10000.;
924  ptB=AliVertexingHFUtils::GetBeautyMotherPt(arrayMC,partDp);
925  }else{
926  pdgCode=-3;
927  }
928  }else{
929  pdgCode=-1;
930  }
931  }
932 
933  Double_t invMass=d->InvMassDplus();
934  Double_t rapid=d->YDplus();
935  fYVsPtNoPid->Fill(ptCand,rapid,invMass);
936  if(passTopolAndPIDCuts) {
937  fYVsPt->Fill(ptCand,rapid,invMass);
938  nOS++;
939  }
940  if(fReadMC && labDp>=0){
941  fYVsPtSigNoPid->Fill(ptCand,rapid, invMass);
942  if(passTopolAndPIDCuts)fYVsPtSig->Fill(ptCand,rapid, invMass);
943  }
944 
945  Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
946  if(isFidAcc){
947 
948  Double_t minPtDau=999.,ptmax=0,maxdca=0,sigvert=0,sumD02=0;
949  Double_t dlen=0,cosp=0,dlenxy=0,cospxy=0, ndlenxy=0, dd0max=0;
951  dlen=d->DecayLength();
952  cosp=d->CosPointingAngle();
953  sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
954  maxdca=-9999.;
955  for(Int_t idau=0;idau<3;idau++) if(d->GetDCA(idau)>maxdca) maxdca=d->GetDCA(idau);
956  sigvert=d->GetSigmaVert();
957  ptmax=0;
958  dlenxy = d->DecayLengthXY();
959  ndlenxy=d->NormalizedDecayLengthXY();
960  cospxy=d->CosPointingAngleXY();
961  for(Int_t i=0; i<3; i++) {
962  if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
963  if(d->PtProng(i)<minPtDau) minPtDau=d->PtProng(i);
964  Double_t diffIP, errdiffIP;
965  d->Getd0MeasMinusExpProng(i,aod->GetMagneticField(),diffIP,errdiffIP);
966  Double_t normdd0= diffIP/errdiffIP;
967  if(i==0) dd0max=normdd0;
968  else if(TMath::Abs(normdd0)>TMath::Abs(dd0max)) dd0max=normdd0;
969  }
970  }
971  Double_t impparXY=d->ImpParXY()*10000.;
972  Double_t resSel=0;
973  if(fRDCutsAnalysis->GetIsSelectedCuts()) resSel=1;
974  if(passTopolAndPIDCuts) resSel=2;
975 
976  //for all THnSParse except for FD
977  Double_t arrayForSparse[kVarForSparse]={invMass,ptCand,impparXY,resSel,minPtDau,sigvert,cosp,cospxy,dlen,dlenxy,ndlenxy,dd0max};
978  //for THnSparse for FD
979  Double_t arrayForSparseFD[kVarForSparseFD]={invMass,ptCand,impparXY,resSel,minPtDau,sigvert,cosp,cospxy,dlen,dlenxy,ndlenxy,dd0max,ptB};
980  Double_t arrayForSparseTrue[kVarForSparseFD]={invMass,ptCand,trueImpParXY,resSel,minPtDau,sigvert,cosp,cospxy,dlen,dlenxy,ndlenxy,dd0max,ptB};
981  Double_t flagOrigin = 0;
982 
983  //Fill histos
984  index=GetHistoIndex(iPtBin);
985  fHistNEvents->Fill(9);
986  nSelectednopid++;
987  fPtVsMassNoPid->Fill(invMass,ptCand);
988  fMassHistNoPid[index]->Fill(invMass);
989  if(fDoImpPar){
990  fHistMassPtImpPar[0]->Fill(arrayForSparse);
991  }
992  if(passTopolAndPIDCuts){
993  fHistNEvents->Fill(10);
994  nSelected++;
995  fPtVsMass->Fill(invMass,ptCand);
996  fMassHist[index]->Fill(invMass);
997  if(d->GetCharge()>0) fMassHistPlus[index]->Fill(invMass);
998  else if(d->GetCharge()<0) fMassHistMinus[index]->Fill(invMass);
999  fPhiEtaCand->Fill(etaD,phiD);
1000  if(TMath::Abs(invMass-1.8696)<0.05) fPhiEtaCandSigReg->Fill(etaD,phiD);
1001  if(fCutsDistr){
1002  fCosPHist[index]->Fill(cosp);
1003  fDLenHist[index]->Fill(dlen);
1004  fSumd02Hist[index]->Fill(sumD02);
1005  fSigVertHist[index]->Fill(sigvert);
1006  fPtMaxHist[index]->Fill(ptmax);
1007  fPtKHist[index]->Fill(d->PtProng(1));
1008  fPtpi1Hist[index]->Fill(d->PtProng(0));
1009  fPtpi2Hist[index]->Fill(d->PtProng(2));
1010  fDCAHist[index]->Fill(maxdca);
1011  fDLxy[index]->Fill(ndlenxy);
1012  fCosxy[index]->Fill(cospxy);
1013  fCorreld0Kd0pi[0]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),d->Getd0Prong(2)*d->Getd0Prong(1));
1014  }
1015  }
1016 
1017  // fill ntuple
1018  if(fFillNtuple==1){
1019 
1020  Float_t tmp[31];
1021 
1022  tmp[0]=pdgCode;
1023  if(isFeeddown) tmp[0]+=5000.;
1024  tmp[1]=d->Px();
1025  tmp[2]=d->Py();
1026  tmp[3]=d->Pz();
1027  tmp[4]=d->Pt();
1028  tmp[5]=d->GetCharge();
1029  // tmp[5]=fRDCutsAnalysis->GetPIDBitMask(d);
1030  tmp[6]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(0));
1031  tmp[7]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(1));
1032  tmp[8]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(2));
1033  tmp[9]=d->PtProng(0);
1034  tmp[10]=d->PtProng(1);
1035  tmp[11]=d->PtProng(2);
1036  tmp[12]=d->PProng(0);
1037  tmp[13]=d->PProng(1);
1038  tmp[14]=d->PProng(2);
1039  tmp[15]=cosp;
1040  tmp[16]=cospxy;
1041  tmp[17]=dlen;
1042  tmp[18]=d->NormalizedDecayLength();
1043  tmp[19]=dlenxy;
1044  tmp[20]=ndlenxy;
1045  tmp[21]=d->InvMassDplus();
1046  tmp[22]=sigvert;
1047  tmp[23]=d->Getd0Prong(0);
1048  tmp[24]=d->Getd0Prong(1);
1049  tmp[25]=d->Getd0Prong(2);
1050  tmp[26]=maxdca;
1051  tmp[27]=ntracks;
1052  tmp[28]=fRDCutsAnalysis->GetCentrality(aod);
1053  tmp[29]=runNumber;
1054  tmp[30]=d->HasBadDaughters();
1055  fNtupleDplus->Fill(tmp);
1056  PostData(4,fNtupleDplus);
1057  }
1058 
1059  if(fFillNtuple==2 && passTopolAndPIDCuts){
1060  Float_t tmp[5];
1061  tmp[0]=pdgCode;
1062  if(isFeeddown) tmp[0]+=5000.;
1063  tmp[1]=d->Pt();
1064  tmp[2]=d->InvMassDplus();
1065  tmp[3]=impparXY;
1066  if(!fReadMC){flagOrigin=0;};
1067  if(fReadMC){
1068  if(isPrimary&&labDp>=0)flagOrigin=1;
1069  if(isFeeddown&&labDp>=0)flagOrigin=2;
1070  if(!(labDp>=0))flagOrigin=3;
1071 
1072  }
1073 
1074  tmp[4]=flagOrigin;
1075  fNtupleDplus->Fill(tmp);
1076  PostData(4,fNtupleDplus);
1077  }
1078 
1079 
1080  if(fReadMC){
1081  if(labDp>=0) {
1082  index=GetSignalHistoIndex(iPtBin);
1083  if(fDoImpPar){
1084  if(isPrimary) fHistMassPtImpPar[1]->Fill(arrayForSparse);
1085  else if(isFeeddown){
1086  fHistMassPtImpPar[2]->Fill(arrayForSparseFD);
1087  fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
1088  }
1089  }
1090  }else{
1091  index=GetBackgroundHistoIndex(iPtBin);
1092  if(fDoImpPar)fHistMassPtImpPar[4]->Fill(arrayForSparse);
1093  }
1094  fMassHistNoPid[index]->Fill(invMass);
1095  if(passTopolAndPIDCuts){
1096  if(fCutsDistr){
1097  Float_t fact=1.;
1098  Float_t factor[3]={1.,1.,1.};
1099  if(fUseStrangeness) fact=GetStrangenessWeights(d,arrayMC,factor);
1100  fCosPHist[index]->Fill(cosp,fact);
1101  fDLenHist[index]->Fill(dlen,fact);
1102  fDLxy[index]->Fill(ndlenxy);
1103  fCosxy[index]->Fill(cospxy);
1104  Float_t sumd02s=d->Getd0Prong(0)*d->Getd0Prong(0)*factor[0]*factor[0]+d->Getd0Prong(1)*d->Getd0Prong(1)*factor[1]*factor[1]+d->Getd0Prong(2)*d->Getd0Prong(2)*factor[2]*factor[2];
1105  fSumd02Hist[index]->Fill(sumd02s);
1106  fSigVertHist[index]->Fill(sigvert,fact);
1107  fPtMaxHist[index]->Fill(ptmax,fact);
1108  fPtKHist[index]->Fill(d->PtProng(1),fact);
1109  fPtpi1Hist[index]->Fill(d->PtProng(0),fact);
1110  fPtpi2Hist[index]->Fill(d->PtProng(2),fact);
1111  fDCAHist[index]->Fill(maxdca,fact);
1112  fCorreld0Kd0pi[1]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),d->Getd0Prong(2)*d->Getd0Prong(1));
1113  }
1114  fMassHist[index]->Fill(invMass);
1115  if(d->GetCharge()>0) fMassHistPlus[index]->Fill(invMass);
1116  else if(d->GetCharge()<0) fMassHistMinus[index]->Fill(invMass);
1117  }
1118  }
1119  }
1120 
1121  if(recVtx)fRDCutsAnalysis->CleanOwnPrimaryVtx(d,aod,origownvtx);
1122 
1123  if(unsetvtx) d->UnsetOwnPrimaryVtx();
1124  }
1125  fCounter->StoreCandidates(aod,nSelectednopid,kTRUE);
1126  fCounter->StoreCandidates(aod,nSelected,kFALSE);
1127  }
1128  delete vHF;
1129  //start LS analysis
1130  if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
1131 
1132  PostData(1,fOutput);
1133  PostData(2,fListCuts);
1134  PostData(3,fCounter);
1135  return;
1136 }
1137 
1138 //________________________________________________________________________
1141 
1142  TString hisname;
1143  Int_t indexLS=0;
1144  Int_t index=0;
1145  Int_t nbins=GetNBinsHistos();
1146  for(Int_t i=0;i<fNPtBins;i++){
1147 
1148  index=GetHistoIndex(i);
1149  indexLS=GetLSHistoIndex(i);
1150 
1151  hisname.Form("hLSPt%d",i);
1152  fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1153  fMassHistLS[indexLS]->Sumw2();
1154 
1155  hisname.Form("hCosPAllPt%dLS",i);
1156  fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1157  fCosPHistLS[index]->Sumw2();
1158  hisname.Form("hDLenAllPt%dLS",i);
1159  fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1160  fDLenHistLS[index]->Sumw2();
1161  hisname.Form("hSumd02AllPt%dLS",i);
1162  fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1163  fSumd02HistLS[index]->Sumw2();
1164  hisname.Form("hSigVertAllPt%dLS",i);
1165  fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1166  fSigVertHistLS[index]->Sumw2();
1167  hisname.Form("hPtMaxAllPt%dLS",i);
1168  fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1169  fPtMaxHistLS[index]->Sumw2();
1170  hisname.Form("hDCAAllPt%dLS",i);
1171  fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1172  fDCAHistLS[index]->Sumw2();
1173 
1174  index=GetSignalHistoIndex(i);
1175  indexLS++;
1176 
1177  hisname.Form("hLSPt%dnw",i);
1178  fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1179  fMassHistLS[indexLS]->Sumw2();
1180 
1181  hisname.Form("hCosPSigPt%dLS",i);
1182  fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1183  fCosPHistLS[index]->Sumw2();
1184  hisname.Form("hDLenSigPt%dLS",i);
1185  fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1186  fDLenHistLS[index]->Sumw2();
1187  hisname.Form("hSumd02SigPt%dLS",i);
1188  fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1189  fSumd02HistLS[index]->Sumw2();
1190  hisname.Form("hSigVertSigPt%dLS",i);
1191  fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1192  fSigVertHistLS[index]->Sumw2();
1193  hisname.Form("hPtMaxSigPt%dLS",i);
1194  fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1195  fPtMaxHistLS[index]->Sumw2();
1196  hisname.Form("hDCASigPt%dLS",i);
1197  fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1198  fDCAHistLS[index]->Sumw2();
1199 
1200  index=GetBackgroundHistoIndex(i);
1201  indexLS++;
1202 
1203  hisname.Form("hLSPt%dLCntrip",i);
1204  fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1205  fMassHistLS[indexLS]->Sumw2();
1206 
1207  hisname.Form("hCosPBkgPt%dLS",i);
1208  fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1209  fCosPHistLS[index]->Sumw2();
1210  hisname.Form("hDLenBkgPt%dLS",i);
1211  fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1212  fDLenHistLS[index]->Sumw2();
1213  hisname.Form("hSumd02BkgPt%dLS",i);
1214  fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1215  fSumd02HistLS[index]->Sumw2();
1216  hisname.Form("hSigVertBkgPt%dLS",i);
1217  fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1218  fSigVertHistLS[index]->Sumw2();
1219  hisname.Form("hPtMaxBkgPt%dLS",i);
1220  fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1221  fPtMaxHistLS[index]->Sumw2();
1222  hisname.Form("hDCABkgPt%dLS",i);
1223  fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1224  fDCAHistLS[index]->Sumw2();
1225 
1226  indexLS++;
1227  hisname.Form("hLSPt%dLCntripsinglecut",i);
1228  fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1229  fMassHistLS[indexLS]->Sumw2();
1230 
1231  indexLS++;
1232  hisname.Form("hLSPt%dspc",i);
1233  fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1234  fMassHistLS[indexLS]->Sumw2();
1235  }
1236 
1237  for(Int_t i=0; i<3*fNPtBins; i++){
1238  fOutput->Add(fCosPHistLS[i]);
1239  fOutput->Add(fDLenHistLS[i]);
1240  fOutput->Add(fSumd02HistLS[i]);
1241  fOutput->Add(fSigVertHistLS[i]);
1242  fOutput->Add(fPtMaxHistLS[i]);
1243  fOutput->Add(fDCAHistLS[i]);
1244  }
1245  for(Int_t i=0; i<5*fNPtBins; i++){
1246  fOutput->Add(fMassHistLS[i]);
1247  }
1248 }
1249 
1250 //________________________________________________________________________
1253 
1254  Int_t nmassbins=GetNBinsHistos();
1255  Double_t maxmult;
1256  if(fSystem==1) maxmult=5000;
1257  else maxmult=200;
1258 
1259  Int_t nptbins=80;
1260  Double_t ptmin=0.;
1261  Double_t ptmax=40.;
1262 
1263  Int_t nptmindaubins=10;
1264  Double_t minptmindau=0.2;
1265  Double_t maxptmindau=1.2;
1266 
1267  Int_t nsigvertbins=25;
1268  Double_t minsigvert=0.010;
1269  Double_t maxsigvert=0.035;
1270 
1271  Int_t nselbins=2;
1272  Double_t minsel=0.5;
1273  Double_t maxsel=2.5;
1274 
1275  Int_t ncospbins=100;
1276  Double_t mincosp=0.90;
1277  Double_t maxcosp=1.;
1278 
1279  Int_t ncospxybins=30;
1280  Double_t mincospxy=0.97;
1281  Double_t maxcospxy=1.;
1282 
1283  Int_t ndeclbins=70;
1284  Double_t mindecl=0.;
1285  Double_t maxdecl=0.7;
1286 
1287  Int_t ndeclxybins=70;
1288  Double_t mindeclxy=0.;
1289  Double_t maxdeclxy=0.7;
1290 
1291  Int_t nnormdlbins=30;
1292  Double_t minnormdl=0.;
1293  Double_t maxnormdl=30.;
1294 
1295  Int_t nd0d0expbins=40;
1296  Double_t mind0d0=-10.;
1297  Double_t maxd0d0=10.;
1298 
1299 
1300  //dimensions for THnSparse which are NOT for BFeed
1301  TString axTit[kVarForSparse]={"M_{K#pi#pi} (GeV/c^{2})","p_{T} (GeV/c)","Imp Par (#mum)","passTopolPID","min. daughter p_{T} (GeV/c)",
1302  "sigmaVertex","cos(#theta_{P})","cos(#theta_{P}^{xy})","decL (cm)","decL XY (cm)","Norm decL XY","Norm max d0-d0exp"};
1303 
1304  Int_t nbins[kVarForSparse]={nmassbins,nptbins,fNImpParBins,nselbins,nptmindaubins,nsigvertbins,ncospbins,ncospxybins,ndeclbins,ndeclxybins,nnormdlbins,nd0d0expbins};
1305  Double_t xmin[kVarForSparse]={fLowmasslimit,ptmin,fLowerImpPar,minsel,minptmindau,minsigvert,mincosp,mincospxy,mindecl,mindeclxy,minnormdl,mind0d0};
1306  Double_t xmax[kVarForSparse]={fUpmasslimit,ptmax,fHigherImpPar,maxsel,maxptmindau,maxsigvert,maxcosp,maxcospxy,maxdecl,maxdeclxy,maxnormdl,maxd0d0};
1307 
1308  //dimensions for THnSparse for BFeed
1309  Int_t nbinsFD[kVarForSparseFD]={nmassbins,nptbins,fNImpParBins,nselbins,nptmindaubins,nsigvertbins,ncospbins,ncospxybins,ndeclbins,ndeclbins,nnormdlbins,nd0d0expbins,84};
1310  Double_t xminFD[kVarForSparseFD]={fLowmasslimit,ptmin,fLowerImpPar,minsel,minptmindau,minsigvert,mincosp,mincospxy,mindecl,mindeclxy,minnormdl,mind0d0,-2};
1311  Double_t xmaxFD[kVarForSparseFD]={fUpmasslimit,ptmax,fHigherImpPar,maxsel,maxptmindau,maxsigvert,maxcosp,maxcospxy,maxdecl,maxdeclxy,maxnormdl,maxd0d0,40};
1312 
1313  //mass, pt, imppar, cosPoinXY, decLXY, norm decLXY (for BFeed also ptB)
1314  //mass, pt, imppar, cosPoinXY, decLXY, norm decLXY (for BFeed also ptB)
1315  fHistMassPtImpPar[0]=new THnSparseF("hMassPtImpParAll",
1316  "Mass vs. pt vs.imppar - All",
1317  kVarForSparse,nbins,xmin,xmax);
1318  fHistMassPtImpPar[1]=new THnSparseF("hMassPtImpParPrompt",
1319  "Mass vs. pt vs.imppar - promptD",
1320  kVarForSparse,nbins,xmin,xmax);
1321  fHistMassPtImpPar[2]=new THnSparseF("hMassPtImpParBfeed",
1322  "Mass vs. pt vs.imppar - DfromB",
1323  kVarForSparseFD,nbinsFD,xminFD,xmaxFD);
1324  fHistMassPtImpPar[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1325  "Mass vs. pt vs.true imppar -DfromB",
1326  kVarForSparseFD,nbinsFD,xminFD,xmaxFD);
1327  fHistMassPtImpPar[4]=new THnSparseF("hMassPtImpParBkg",
1328  "Mass vs. pt vs.imppar - backgr.",
1329  kVarForSparse,nbins,xmin,xmax);
1330 
1331 
1332  for(Int_t i=0; i<5; i++){
1333  for(Int_t iax=0; iax<kVarForSparse; iax++) fHistMassPtImpPar[i]->GetAxis(iax)->SetTitle(axTit[iax].Data());
1334  if(i == 2 || i == 3) fHistMassPtImpPar[i]->GetAxis(kVarForSparseFD-1)->SetTitle("p_{T}^{B} (GeV/c)");
1335  fOutput->Add(fHistMassPtImpPar[i]);
1336  }
1337 }
1338 
1339 //________________________________________________________________________
1342 
1343  const Int_t nVarPrompt = 2;
1344  const Int_t nVarFD = 3;
1345 
1346  Int_t nbinsPrompt[nVarPrompt]={200,100};
1347  Int_t nbinsFD[nVarFD]={200,100,200};
1348 
1349  Double_t xminPrompt[nVarPrompt] = {0.,-1.};
1350  Double_t xmaxPrompt[nVarPrompt] = {40.,1.};
1351 
1352  Double_t xminFD[nVarFD] = {0.,-1.,0.};
1353  Double_t xmaxFD[nVarFD] = {40.,1.,40.};
1354 
1355  //pt, y
1356  fMCAccPrompt = new THnSparseF("hMCAccPrompt","kStepMCAcceptance pt vs. y - promptD",nVarPrompt,nbinsPrompt,xminPrompt,xmaxPrompt);
1357  fMCAccPrompt->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
1358  fMCAccPrompt->GetAxis(1)->SetTitle("y");
1359 
1360  //pt,y,ptB
1361  fMCAccBFeed = new THnSparseF("hMCAccBFeed","kStepMCAcceptance pt vs. y vs. ptB - DfromB",nVarFD,nbinsFD,xminFD,xmaxFD);
1362  fMCAccBFeed->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
1363  fMCAccBFeed->GetAxis(1)->SetTitle("y");
1364  fMCAccBFeed->GetAxis(2)->SetTitle("p_{T}^{B} (GeV/c)");
1365 
1366  fOutput->Add(fMCAccPrompt);
1367  fOutput->Add(fMCAccBFeed);
1368 }
1369 
1370 //________________________________________________________________________
1371 void AliAnalysisTaskSEDplus::FillMCAcceptanceHistos(TClonesArray *arrayMC, AliAODMCHeader *mcHeader){
1373 
1374  const Int_t nProng = 3;
1375 
1376  Double_t zMCVertex = mcHeader->GetVtxZ(); //vertex MC
1377 
1378  for(Int_t iPart=0; iPart<arrayMC->GetEntriesFast(); iPart++){
1379  AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(arrayMC->At(iPart));
1380  if (TMath::Abs(mcPart->GetPdgCode()) == 411){
1381 
1382  Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,mcPart,fUseQuarkTagInKine);//Prompt = 4, FeedDown = 5
1383 
1384  Int_t deca = 0;
1385  Bool_t isGoodDecay=kFALSE;
1386  Int_t labDau[4]={-1,-1,-1,-1};
1387  Bool_t isInAcc = kFALSE;
1388  Bool_t isFidAcc = kFALSE;
1389 
1390  deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,mcPart,labDau);
1391  if(deca > 0) isGoodDecay=kTRUE;
1392 
1393  if(labDau[0]==-1){
1394  continue; //protection against unfilled array of labels
1395  }
1396 
1397  isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(mcPart->Pt(),mcPart->Y());
1398  isInAcc=CheckAcc(arrayMC,nProng,labDau);
1399 
1400  if(isGoodDecay && TMath::Abs(zMCVertex) < fRDCutsAnalysis->GetMaxVtxZ() && isFidAcc && isInAcc) {
1401  //for prompt
1402  if(orig == 4){
1403  //fill histo for prompt
1404  Double_t arrayMCprompt[2] = {mcPart->Pt(),mcPart->Y()};
1405  fMCAccPrompt->Fill(arrayMCprompt);
1406  }
1407  //for FD
1408  else if(orig == 5){
1409  Double_t ptB = AliVertexingHFUtils::GetBeautyMotherPt(arrayMC,mcPart);
1410  //fill histo for FD
1411  Double_t arrayMCFD[3] = {mcPart->Pt(),mcPart->Y(),ptB};
1412  fMCAccBFeed->Fill(arrayMCFD);
1413  }
1414  else
1415  continue;
1416  }
1417  }
1418  }
1419 }
1420 
1421 //________________________________________________________________________
1422 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
1423 {
1425  //
1426  if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
1427 
1428  fOutput = dynamic_cast<TList*> (GetOutputData(1));
1429  if (!fOutput) {
1430  printf("ERROR: fOutput not available\n");
1431  return;
1432  }
1433 
1434  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
1435  if(fHistNEvents){
1436  printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
1437  }else{
1438  printf("ERROR: fHistNEvents not available\n");
1439  return;
1440  }
1441 
1442  return;
1443 }
1444 //_________________________________________________________________________________________________
1445 Float_t AliAnalysisTaskSEDplus::GetTrueImpactParameter(const AliAODMCHeader *mcHeader, TClonesArray* arrayMC, const AliAODMCParticle *partDp) const {
1447 
1448  Double_t vtxTrue[3];
1449  mcHeader->GetVertex(vtxTrue);
1450  Double_t origD[3];
1451  partDp->XvYvZv(origD);
1452  Short_t charge=partDp->Charge();
1453  Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
1454  for(Int_t iDau=0; iDau<3; iDau++){
1455  pXdauTrue[iDau]=0.;
1456  pYdauTrue[iDau]=0.;
1457  pZdauTrue[iDau]=0.;
1458  }
1459 
1460  Int_t nDau=partDp->GetNDaughters();
1461  Int_t labelFirstDau = partDp->GetDaughter(0);
1462  if(nDau==3){
1463  for(Int_t iDau=0; iDau<3; iDau++){
1464  Int_t ind = labelFirstDau+iDau;
1465  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1466  if(!part){
1467  AliError("Daughter particle not found in MC array");
1468  return 99999.;
1469  }
1470  pXdauTrue[iDau]=part->Px();
1471  pYdauTrue[iDau]=part->Py();
1472  pZdauTrue[iDau]=part->Pz();
1473  }
1474  }else if(nDau==2){
1475  Int_t theDau=0;
1476  for(Int_t iDau=0; iDau<2; iDau++){
1477  Int_t ind = labelFirstDau+iDau;
1478  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1479  if(!part){
1480  AliError("Daughter particle not found in MC array");
1481  return 99999.;
1482  }
1483  Int_t pdgCode=TMath::Abs(part->GetPdgCode());
1484  if(pdgCode==211 || pdgCode==321){
1485  pXdauTrue[theDau]=part->Px();
1486  pYdauTrue[theDau]=part->Py();
1487  pZdauTrue[theDau]=part->Pz();
1488  ++theDau;
1489  }else{
1490  Int_t nDauRes=part->GetNDaughters();
1491  if(nDauRes==2){
1492  Int_t labelFirstDauRes = part->GetDaughter(0);
1493  for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
1494  Int_t indDR = labelFirstDauRes+iDauRes;
1495  AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
1496  if(!partDR){
1497  AliError("Daughter particle not found in MC array");
1498  return 99999.;
1499  }
1500 
1501  Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
1502  if(pdgCodeDR==211 || pdgCodeDR==321){
1503  pXdauTrue[theDau]=partDR->Px();
1504  pYdauTrue[theDau]=partDR->Py();
1505  pZdauTrue[theDau]=partDR->Pz();
1506  ++theDau;
1507  }
1508  }
1509  }
1510  }
1511  }
1512  if(theDau!=3){
1513  AliError("Wrong number of decay prongs");
1514  return 99999.;
1515  }
1516  }
1517 
1518  Double_t d0dummy[3]={0.,0.,0.};
1519  AliAODRecoDecayHF aodDplusMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
1520  return aodDplusMC.ImpParXY();
1521 
1522 }
1523 //_________________________________________________________________________________________________
1524 Float_t AliAnalysisTaskSEDplus::GetStrangenessWeights(const AliAODRecoDecayHF3Prong* d, TClonesArray* arrayMC, Float_t factor[3]) const {
1526 
1527  for(Int_t iprong=0;iprong<3;iprong++){
1528  factor[iprong]=1;
1529  AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
1530  Int_t labd= trad->GetLabel();
1531  if(labd>=0){
1532  AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
1533  if(dau){
1534  Int_t labm = dau->GetMother();
1535  if(labm>=0){
1536  AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
1537  if(mot){
1538  if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
1539  if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
1540  else factor[iprong]=1./.6;
1541  // fNentries->Fill(12);
1542  }
1543  if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
1544  factor[iprong]=1./0.25;
1545  // fNentries->Fill(13);
1546  }//if 3122
1547  }//if(mot)
1548  }//if labm>0
1549  }//if(dau)
1550  }//if labd>=0
1551  }//prong loop
1552 
1553  Float_t fact=1.;
1554  for(Int_t k=0;k<3;k++)fact=fact*factor[k];
1555  return fact;
1556 
1557 }
1558 
1559 //_________________________________________________________________
1560 Bool_t AliAnalysisTaskSEDplus::CheckAcc(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
1562  for (Int_t iProng = 0; iProng<nProng; iProng++){
1563  AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
1564  if(!mcPartDaughter) return kFALSE;
1565  Double_t eta = mcPartDaughter->Eta();
1566  Double_t pt = mcPartDaughter->Pt();
1567  if (TMath::Abs(eta) > 0.9 || pt < 0.1) return kFALSE;
1568  }
1569  return kTRUE;
1570 }
TH1F * fCosxy[3 *kMaxPtBins]
!hist. for Cosxy (topol+PID)
Int_t charge
THnSparseF * fMCAccPrompt
!histo for StepMCAcc for Dplus prompt (pt,y,ptB)
TH2F * fPhiEtaCand
! hist. with eta/phi distribution of candidates
Double_t NormalizedDecayLengthXY() const
Double_t NormalizedDecayLength() const
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
Bool_t fUseStrangeness
flag for access to MC
TH1F * fSigVertHist[3 *kMaxPtBins]
!hist. for sigVert (topol+PID)
static Int_t CheckDplusDecay(AliStack *stack, Int_t label, Int_t *arrayDauLab)
AliRDHFCutsDplustoKpipi * fRDCutsAnalysis
list of cuts
virtual void UserCreateOutputObjects()
Implementation of interface methods.
void StoreCandidates(AliVEvent *, Int_t nCand=0, Bool_t flagFilter=kTRUE)
THnSparseF * fHistMassPtImpPar[5]
! histograms for impact parameter and cut variation study
TH2F * fYVsPtSigNoPid
! hist. of Y vs. Pt (MC, only sig, w/o PID)
void Getd0MeasMinusExpProng(Int_t ip, Double_t magf, Double_t &d0diff, Double_t &errd0diff) const
TH1F * fPtpi1Hist[3 *kMaxPtBins]
!hist. for PtPi1 (topol+PID)
Int_t GetIsSelectedCuts() const
Definition: AliRDHFCuts.h:334
Int_t IsEventSelectedInCentrality(AliVEvent *event)
Bool_t HasSelectionBit(Int_t i) const
TH3F * fYVsPt
! hist. of Y vs. Pt vs. Mass (topol+PID cuts)
virtual void Terminate(Option_t *option)
TH1F * fDLxy[3 *kMaxPtBins]
!hist. for DLxy (topol+PID)
Int_t fEtaSelection
flag to do LS analysis
Float_t fLowerImpPar
nunber of bins in impact parameter histos
TH2F * fHistCentrality[3]
!hist. for cent distr (all,sel ev, )
TH1F * fPtMaxHistLS[3 *kMaxPtBins]
!hist. for LS cuts variable 5 (topol+PID)
Double_t ImpParXY() const
TH1F * fSPDMult
! hist. of spd mult
Bool_t fReadMC
flag for filling ntuple 0 no NTuple 1 big Ntuple 2 small NTuple
Int_t GetLSHistoIndex(Int_t iPtBin) const
Int_t GetSignalHistoIndex(Int_t iPtBin) const
UInt_t GetPIDTrackTPCTOFBitMap(AliAODTrack *track) const
Int_t GetWhyRejection() const
Definition: AliRDHFCuts.h:294
Double_t CosPointingAngleXY() const
TH1F * fMassHistNoPid[3 *kMaxPtBins]
!hist. for inv mass (w/o PID)
TList * fListCuts
width of one bin in output histos
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
TList * fOutput
! list send on output slot 0
virtual Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y) const
Bool_t fCutsDistr
flag to use bitmask
static Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark=kTRUE)
Double_t GetMaxVtxZ() const
Definition: AliRDHFCuts.h:241
void LSAnalysis(TClonesArray *arrayOppositeSign, TClonesArray *arrayLikeSign, AliAODEvent *aod, AliAODVertex *vtx1, Int_t nDplusOS)
Int_t GetBackgroundHistoIndex(Int_t iPtBin) const
TH1F * fCosPHist[3 *kMaxPtBins]
!hist. for PointingAngle (topol+PID)
Bool_t HasBadDaughters() const
TH2F * fYVsPtSig
! hist. of Y vs. Pt (MC, only sig, topol+PID cuts)
Float_t fHigherImpPar
lower limit in impact parameter (um)
Class for cuts on AOD reconstructed D+->Kpipi.
TH1F * fMassHistLS[5 *kMaxPtBins]
!hist. for LS inv mass (topol+PID)
static Int_t GetNumberOfTrackletsInEtaRange(AliAODEvent *ev, Double_t mineta, Double_t maxeta)
Int_t fFillNtuple
limits for the Pt bins
TH1F * fSumd02Hist[3 *kMaxPtBins]
!hist. for sum d02 (topol+PID)
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
const Double_t ptmax
void FillMCAcceptanceHistos(TClonesArray *arrayMC, AliAODMCHeader *mcHeader)
TH1F * fMassHistPlus[3 *kMaxPtBins]
!hist. for D+ inv mass (topol+PID cuts)
TH2F * fPtVsMass
! hist. of pt vs. mass (topol+PID cuts)
AliAODVertex * GetOwnPrimaryVtx() const
Double_t GetSigmaVert(const AliAODEvent *aod=0x0)
Int_t fDoLS
higher limit in impact parameter (um)
TH1F * fPtpi2Hist[3 *kMaxPtBins]
!hist. for PtPi2 (topol+PID)
TH2F * fPhiEtaCandSigReg
! hist. eta/phi of candidates in D+ mass region
const Double_t ptmin
Float_t GetCentrality(AliAODEvent *aodEvent)
Definition: AliRDHFCuts.h:242
Float_t GetTrueImpactParameter(const AliAODMCHeader *mcHeader, TClonesArray *arrayMC, const AliAODMCParticle *partDp) const
Int_t fNImpParBins
flag for quark/hadron level identification of prompt and feeddown
TH1F * fDLenHistLS[3 *kMaxPtBins]
!hist. for LS cuts variable 2 (topol+PID)
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
TH3F * fYVsPtNoPid
! hist. of Y vs. Pt vs. Mass(w/o PID)
Int_t GetHistoIndex(Int_t iPtBin) const
THnSparseF * fMCAccBFeed
!histo for StepMCAcc for Dplus FD (pt,y,ptB)
TH1F * fMassHist[3 *kMaxPtBins]
!hist. for inv mass (topol+PID cuts)
void SetOwnPrimaryVtx(const AliAODVertex *vtx)
TH2F * fPtVsMassNoPid
! hist. of pt vs. mass (w/o PID)
Double_t DecayLengthXY() const
static Double_t GetBeautyMotherPt(TClonesArray *arrayMC, AliAODMCParticle *mcPart)
Bool_t CheckAcc(TClonesArray *arrayMC, Int_t nProng, Int_t *labDau)
Bool_t GetIsPrimaryWithoutDaughters() const
Definition: AliRDHFCuts.h:262
TH1F * fSumd02HistLS[3 *kMaxPtBins]
!hist. for LS cuts variable 3 (topol+PID)
Bool_t IsEventSelected(AliVEvent *event)
void StoreEvent(AliVEvent *, AliRDHFCuts *, Bool_t mc=kFALSE, Int_t multiplicity=-9999, Double_t spherocity=-99.)
TH1F * fHistNEvents
!hist. for No. of events
TH1F * fCosPHistLS[3 *kMaxPtBins]
!hist. for LS cuts variable 1 (topol+PID)
TH1F * fPtKHist[3 *kMaxPtBins]
!hist. for PtK (topol+PID)
Float_t fBinWidth
Number of Pt Bins.
void SetMassLimits(Float_t range)
void CleanOwnPrimaryVtx(AliAODRecoDecayHF *d, AliAODEvent *aod, AliAODVertex *origownvtx) const
Int_t fSystem
eta region to accept D+ 0=all, -1 = negative, 1 = positive
TH2F * fCorreld0Kd0pi[3]
!hist. for d0k*d0pi vs. d0k*d0pi (topol+PID)
virtual void UserExec(Option_t *option)
TH1F * fDLenHist[3 *kMaxPtBins]
!hist. for Dec Length (topol+PID)
Bool_t fDoImpPar
flag to activate cuts distr histos
AliNormalizationCounter * fCounter
Cuts for Analysis.
TH1F * fSigVertHistLS[3 *kMaxPtBins]
!hist. for LS cuts variable 4 (topol+PID)
Bool_t fUseBit
flag to enhance strangeness in MC to fit to data
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:233
Double_t fArrayBinLimits[kMaxPtBins+1]
const Int_t nbins
Double_t CosPointingAngle() const
Float_t GetStrangenessWeights(const AliAODRecoDecayHF3Prong *d, TClonesArray *arrayMC, Float_t factor[3]) const
TH1F * fDCAHist[3 *kMaxPtBins]
!hist. for DCA (topol+PID)
Float_t fLowmasslimit
upper inv mass limit for histos
TH1F * fMassHistMinus[3 *kMaxPtBins]
!hist. for D- inv mass (topol+PID cuts)
Int_t GetUseCentrality() const
Definition: AliRDHFCuts.h:264
Bool_t RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d, AliAODEvent *aod) const
Int_t PtBin(Double_t pt) const
Double_t DecayLength() const
TNtuple * fNtupleDplus
! output ntuple
Bool_t fStepMCAcc
flag to activate impact paramter histos
Int_t fNPtBins
lower inv mass limit for histos
Bool_t fUseQuarkTagInKine
flag to activate histos for StepMCAcc
Int_t nptbins
TH1F * fDCAHistLS[3 *kMaxPtBins]
!hist. for LS cuts variable 6 (topol+PID)
TH1F * fPtMaxHist[3 *kMaxPtBins]
!hist. for Pt Max (topol+PID)