AliPhysics  29d4213 (29d4213)
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Macros
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 ",11,-0.5,10.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 
477  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
478  fHistNEvents->Sumw2();
479  fHistNEvents->SetMinimum(0);
480  fOutput->Add(fHistNEvents);
481 
482  TString hisname;
483  Int_t index=0;
484  Int_t nbins=GetNBinsHistos();
485  fHistCentrality[0]=new TH2F("hCentrMult","centrality",100,0.5,30000.5,40,0.,100.);
486  fHistCentrality[1]=new TH2F("hCentrMult(selectedCent)","centrality(selectedCent)",100,0.5,30000.5,40,0.,100.);
487  fHistCentrality[2]=new TH2F("hCentrMult(OutofCent)","centrality(OutofCent)",100,0.5,30000.5,40,0.,100.);
488  for(Int_t i=0;i<3;i++){
489  fHistCentrality[i]->Sumw2();
490  fOutput->Add(fHistCentrality[i]);
491  }
492  for(Int_t i=0;i<fNPtBins;i++){
493 
494  index=GetHistoIndex(i);
495 
496  hisname.Form("hMassNoPidPt%d",i);
497  fMassHistNoPid[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
498  fMassHistNoPid[index]->Sumw2();
499  hisname.Form("hCosPAllPt%d",i);
500  fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
501  fCosPHist[index]->Sumw2();
502  hisname.Form("hDLenAllPt%d",i);
503  fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
504  fDLenHist[index]->Sumw2();
505  hisname.Form("hSumd02AllPt%d",i);
506  fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
507  fSumd02Hist[index]->Sumw2();
508  hisname.Form("hSigVertAllPt%d",i);
509  fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
510  fSigVertHist[index]->Sumw2();
511  hisname.Form("hPtMaxAllPt%d",i);
512  fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
513  fPtMaxHist[index]->Sumw2();
514  hisname.Form("hPtKPt%d",i);
515  fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
516  fPtKHist[index]->Sumw2();
517  hisname.Form("hPtpi1Pt%d",i);
518  fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
519  fPtpi1Hist[index]->Sumw2();
520  hisname.Form("hPtpi2Pt%d",i);
521  fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
522  fPtpi2Hist[index]->Sumw2();
523  hisname.Form("hDCAAllPt%d",i);
524  fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
525  fDCAHist[index]->Sumw2();
526 
527  hisname.Form("hDLxyPt%d",i);
528  fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
529  fDLxy[index]->Sumw2();
530  hisname.Form("hCosxyPt%d",i);
531  fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
532  fCosxy[index]->Sumw2();
533 
534  hisname.Form("hMassPt%dTC",i);
535  fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
536  fMassHist[index]->Sumw2();
537  hisname.Form("hMassPt%dTCPlus",i);
538  fMassHistPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
539  fMassHistPlus[index]->Sumw2();
540  hisname.Form("hMassPt%dTCMinus",i);
541  fMassHistMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
542  fMassHistMinus[index]->Sumw2();
543 
544 
545 
546  index=GetSignalHistoIndex(i);
547  hisname.Form("hSigNoPidPt%d",i);
548  fMassHistNoPid[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
549  fMassHistNoPid[index]->Sumw2();
550  hisname.Form("hCosPSigPt%d",i);
551  fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
552  fCosPHist[index]->Sumw2();
553  hisname.Form("hDLenSigPt%d",i);
554  fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
555  fDLenHist[index]->Sumw2();
556  hisname.Form("hSumd02SigPt%d",i);
557  fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
558  fSumd02Hist[index]->Sumw2();
559  hisname.Form("hSigVertSigPt%d",i);
560  fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
561  fSigVertHist[index]->Sumw2();
562  hisname.Form("hPtMaxSigPt%d",i);
563  fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
564  fPtMaxHist[index]->Sumw2();
565  hisname.Form("hPtKSigPt%d",i);
566  fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
567  fPtKHist[index]->Sumw2();
568  hisname.Form("hPtpi1SigPt%d",i);
569  fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
570  fPtpi1Hist[index]->Sumw2();
571  hisname.Form("hPtpi2SigPt%d",i);
572  fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
573  fPtpi2Hist[index]->Sumw2();
574 
575  hisname.Form("hDCASigPt%d",i);
576  fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
577  fDCAHist[index]->Sumw2();
578 
579  hisname.Form("hDLxySigPt%d",i);
580  fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
581  fDLxy[index]->Sumw2();
582  hisname.Form("hCosxySigPt%d",i);
583  fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
584  fCosxy[index]->Sumw2();
585  hisname.Form("hSigPt%dTC",i);
586  fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
587  fMassHist[index]->Sumw2();
588  hisname.Form("hSigPt%dTCPlus",i);
589  fMassHistPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
590  fMassHistPlus[index]->Sumw2();
591  hisname.Form("hSigPt%dTCMinus",i);
592  fMassHistMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
593  fMassHistMinus[index]->Sumw2();
594 
595 
596  index=GetBackgroundHistoIndex(i);
597  hisname.Form("hBkgNoPidPt%d",i);
598  fMassHistNoPid[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
599  fMassHistNoPid[index]->Sumw2();
600  hisname.Form("hCosPBkgPt%d",i);
601  fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
602  fCosPHist[index]->Sumw2();
603  hisname.Form("hDLenBkgPt%d",i);
604  fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
605  fDLenHist[index]->Sumw2();
606  hisname.Form("hSumd02BkgPt%d",i);
607  fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
608  fSumd02Hist[index]->Sumw2();
609  hisname.Form("hSigVertBkgPt%d",i);
610  fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
611  fSigVertHist[index]->Sumw2();
612  hisname.Form("hPtMaxBkgPt%d",i);
613  fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
614  fPtMaxHist[index]->Sumw2();
615  hisname.Form("hPtKBkgPt%d",i);
616  fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
617  fPtKHist[index]->Sumw2();
618  hisname.Form("hPtpi1BkgPt%d",i);
619  fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
620  fPtpi1Hist[index]->Sumw2();
621  hisname.Form("hPtpi2BkgPt%d",i);
622  fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
623  fPtpi2Hist[index]->Sumw2();
624  hisname.Form("hDCABkgPt%d",i);
625  fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
626  fDCAHist[index]->Sumw2();
627 
628  hisname.Form("hDLxyBkgPt%d",i);
629  fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
630  fDLxy[index]->Sumw2();
631  hisname.Form("hCosxyBkgPt%d",i);
632  fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
633  fCosxy[index]->Sumw2();
634 
635 
636  hisname.Form("hBkgPt%dTC",i);
637  fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
638  fMassHist[index]->Sumw2();
639  hisname.Form("hBkgPt%dTCPlus",i);
640  fMassHistPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
641  fMassHistPlus[index]->Sumw2();
642  hisname.Form("hBkgPt%dTCMinus",i);
643  fMassHistMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
644  fMassHistMinus[index]->Sumw2();
645  }
646 
647 
648  for(Int_t i=0; i<3*fNPtBins; i++){
649  fOutput->Add(fMassHistNoPid[i]);
650  if(fCutsDistr){
651  fOutput->Add(fCosPHist[i]);
652  fOutput->Add(fDLenHist[i]);
653  fOutput->Add(fSumd02Hist[i]);
654  fOutput->Add(fSigVertHist[i]);
655  fOutput->Add(fPtMaxHist[i]);
656  fOutput->Add(fPtKHist[i]);
657  fOutput->Add(fPtpi1Hist[i]);
658  fOutput->Add(fPtpi2Hist[i]);
659  fOutput->Add(fDCAHist[i]);
660  fOutput->Add(fDLxy[i]);
661  fOutput->Add(fCosxy[i]);
662  }
663  fOutput->Add(fMassHist[i]);
664  fOutput->Add(fMassHistPlus[i]);
665  fOutput->Add(fMassHistMinus[i]);
666 
667  }
668 
669  fCorreld0Kd0pi[0]=new TH2F("hCorreld0Kd0piAll","",100,-0.02,0.02,100,-0.02,0.02);
670  fCorreld0Kd0pi[1]=new TH2F("hCorreld0Kd0piSig","",100,-0.02,0.02,100,-0.02,0.02);
671  fCorreld0Kd0pi[2]=new TH2F("hCorreld0Kd0piBkg","",100,-0.02,0.02,100,-0.02,0.02);
672  if(fCutsDistr){
673  for(Int_t i=0; i<3; i++){
674  fCorreld0Kd0pi[i]->Sumw2();
675  fOutput->Add(fCorreld0Kd0pi[i]);
676  }
677  }
678 
679 
680  fPtVsMassNoPid=new TH2F("hPtVsMassNoPid","PtVsMass (no PID)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
681  fPtVsMass=new TH2F("hPtVsMass","PtVsMass",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
682  fYVsPtNoPid=new TH3F("hYVsPtNoPid","YvsPt (no PID)",40,0.,20.,80,-2.,2.,nbins,fLowmasslimit,fUpmasslimit);
683  fYVsPt=new TH3F("hYVsPt","YvsPt",40,0.,20.,80,-2.,2.,nbins,fLowmasslimit,fUpmasslimit);
684  fYVsPtSigNoPid=new TH2F("hYVsPtSigNoPid","YvsPt (MC, only sig., no PID)",40,0.,20.,80,-2.,2.);
685  fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only Sig)",40,0.,20.,80,-2.,2.);
686  fPhiEtaCand=new TH2F("hPhiEtaCand","phi vs. eta candidates",20,-1.,1.,50,0.,2*TMath::Pi());
687  fPhiEtaCandSigReg=new TH2F("hPhiEtaCandSigReg","phi vs. eta candidates",20,-1.,1.,50,0.,2*TMath::Pi());
688 
689  Double_t maxmult;
690  if(fSystem==1) maxmult=5000;
691  else maxmult=200;
692  fSPDMult = new TH1F("hSPDMult", "Tracklets multiplicity; Tracklets ; Entries",200,0.,maxmult);
693  fOutput->Add(fPtVsMassNoPid);
694  fOutput->Add(fPtVsMass);
695  fOutput->Add(fYVsPtNoPid);
696  fOutput->Add(fYVsPt);
697  fOutput->Add(fYVsPtSigNoPid);
698  fOutput->Add(fYVsPtSig);
699  fOutput->Add(fPhiEtaCand);
701  fOutput->Add(fSPDMult);
702 
703 
704  //Counter for Normalization
705  TString normName="NormalizationCounter";
706  AliAnalysisDataContainer *cont = GetOutputSlot(3)->GetContainer();
707  if(cont)normName=(TString)cont->GetName();
708  fCounter = new AliNormalizationCounter(normName.Data());
709  fCounter->Init();
710 
714 
715  if(fFillNtuple==1){
716  OpenFile(4); // 4 is the slot number of the ntuple
717 
718 
719  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");
720 
721  }
722 
723  if(fFillNtuple==2){
724  OpenFile(4); // 4 is the slot number of the ntuple
725 
726 
727  fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Pt:InvMass:d0:origin");
728 
729  }
730 
731  return;
732 }
733 
734 //________________________________________________________________________
735 void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
736 {
739 
740  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
741 
742  TClonesArray *array3Prong = 0;
743  TClonesArray *arrayLikeSign =0;
744  if(!aod && AODEvent() && IsStandardAOD()) {
745  // In case there is an AOD handler writing a standard AOD, use the AOD
746  // event in memory rather than the input (ESD) event.
747  aod = dynamic_cast<AliAODEvent*> (AODEvent());
748  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
749  // have to taken from the AOD event hold by the AliAODExtension
750  AliAODHandler* aodHandler = (AliAODHandler*)
751  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
752  if(aodHandler->GetExtensions()) {
753  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
754  AliAODEvent *aodFromExt = ext->GetAOD();
755  array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
756  arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
757  }
758  } else if(aod) {
759  array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
760  arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
761  }
762 
763  if(!aod || !array3Prong) {
764  printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
765  return;
766  }
767  if(!arrayLikeSign && fDoLS) {
768  printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
769  return;
770  }
771 
772  // fix for temporary bug in ESDfilter
773  // the AODs with null vertex pointer didn't pass the PhysSel
774  if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
775 
776  //Store the event in AliNormalizationCounter->To disable for Pb-Pb? Add a flag to disable it?
778 
779  fHistNEvents->Fill(0); // count event
780  Int_t runNumber=aod->GetRunNumber();
781 
782  //Event selection
783  Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
784  Float_t ntracks=aod->GetNumberOfTracks();
785  Float_t evCentr=0;
787  fHistCentrality[0]->Fill(ntracks,evCentr);
788  if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(2);
789  if(fRDCutsAnalysis->GetWhyRejection()==1) fHistNEvents->Fill(3);
790  if(fRDCutsAnalysis->GetWhyRejection()==2){fHistNEvents->Fill(4);fHistCentrality[2]->Fill(ntracks,evCentr);}
791  if(fRDCutsAnalysis->GetWhyRejection()==6)fHistNEvents->Fill(5);
792  if(fRDCutsAnalysis->GetWhyRejection()==7)fHistNEvents->Fill(6);
793 
794  TClonesArray *arrayMC=0;
795  AliAODMCHeader *mcHeader=0;
796  // load MC particles
797  if(fReadMC){
798  arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
799  if(!arrayMC) {
800  printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
801  return;
802  }
803 
804  // load MC header
805  mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
806  if(!mcHeader) {
807  printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
808  return;
809  }
810  }
811  if(fReadMC && fStepMCAcc){
812  if(aod->GetTriggerMask()==0 &&
813  (runNumber>=195344 && runNumber<=195677)){
814  // protection for events with empty trigger mask in p-Pb
815  return;
816  }
818  // events not passing the centrality selection can be removed immediately.
819 
820  FillMCAcceptanceHistos(arrayMC, mcHeader);
821  }
822  // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
823  //TString trigclass=aod->GetFiredTriggerClasses();
824  // Post the data already here
825  PostData(1,fOutput);
826  if(!isEvSel)return;
827  Int_t tracklets=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
828  // printf("ntracklet===%d\n",tracklets);
829  fSPDMult->Fill(tracklets);
830 
831  fHistCentrality[1]->Fill(ntracks,evCentr);
832  fHistNEvents->Fill(1);
833 
834  // AOD primary vertex
835  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
836  // vtx1->Print();
837  // TString primTitle = vtx1->GetTitle();
838  //if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0)fHistNEvents->Fill(2);
839 
840 
841  Int_t n3Prong = array3Prong->GetEntriesFast();
842  // printf("Number of D+->Kpipi: %d and of tracks: %d\n",n3Prong,aod->GetNumberOfTracks());
843 
844  Int_t nOS=0;
845  Int_t index;
846  Int_t pdgDgDplustoKpipi[3]={321,211,211};
847 
848  if(fDoLS>1){//Normalizations for LS
849  for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
850  AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
853  }
854  }
855  }else{//Standard analysis
856  // 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
857  //Double_t *cutsDplus = new (Double_t*)fRDCuts->GetCuts();
858  Int_t nSelectednopid=0,nSelected=0;
859  for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
860  AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
861  fHistNEvents->Fill(7);
863  fHistNEvents->Fill(8);
864  continue;
865  }
866 
867  Int_t passTopolAndPIDCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
868 
869  if(!fRDCutsAnalysis->GetIsSelectedCuts()) continue;
870 
871  Double_t etaD=d->Eta();
872  Double_t phiD=d->Phi();
873  if(fEtaSelection!=0){
874  if(fEtaSelection==1 && etaD<0) continue;
875  if(fEtaSelection==-1 && etaD>0) continue;
876  }
877 
878  Bool_t unsetvtx=kFALSE;
879  if(!d->GetOwnPrimaryVtx()){
880  d->SetOwnPrimaryVtx(vtx1);
881  unsetvtx=kTRUE;
882  }
883 
884  Double_t ptCand = d->Pt();
885  Int_t iPtBin = fRDCutsAnalysis->PtBin(ptCand);
886 
887  Bool_t recVtx=kFALSE;
888  AliAODVertex *origownvtx=0x0;
890  if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
891  if(fRDCutsAnalysis->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
892  else fRDCutsAnalysis->CleanOwnPrimaryVtx(d,aod,origownvtx);
893  }
894 
895  Int_t labDp=-1;
896  Bool_t isPrimary=kFALSE;
897  Bool_t isFeeddown=kFALSE;
898  Float_t pdgCode=-2;
899  Float_t trueImpParXY=0.;
900  Double_t ptB=-1.5;
901  if(fReadMC){
902  labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
903  if(labDp>=0){
904  AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
905  Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,partDp,fUseQuarkTagInKine);//Prompt = 4, FeedDown = 5
906  pdgCode=TMath::Abs(partDp->GetPdgCode());
907  if(orig==4){
908  isPrimary=kTRUE;
909  isFeeddown=kFALSE;
910  }else if(orig==5){
911  isPrimary=kFALSE;
912  isFeeddown=kTRUE;
913  trueImpParXY=GetTrueImpactParameter(mcHeader,arrayMC,partDp)*10000.;
914  ptB=AliVertexingHFUtils::GetBeautyMotherPt(arrayMC,partDp);
915  }else{
916  pdgCode=-3;
917  }
918  }else{
919  pdgCode=-1;
920  }
921  }
922 
923  Double_t invMass=d->InvMassDplus();
924  Double_t rapid=d->YDplus();
925  fYVsPtNoPid->Fill(ptCand,rapid,invMass);
926  if(passTopolAndPIDCuts) {
927  fYVsPt->Fill(ptCand,rapid,invMass);
928  nOS++;
929  }
930  if(fReadMC && labDp>=0){
931  fYVsPtSigNoPid->Fill(ptCand,rapid, invMass);
932  if(passTopolAndPIDCuts)fYVsPtSig->Fill(ptCand,rapid, invMass);
933  }
934 
935  Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
936  if(isFidAcc){
937 
938  Double_t minPtDau=999.,ptmax=0,maxdca=0,sigvert=0,sumD02=0;
939  Double_t dlen=0,cosp=0,dlenxy=0,cospxy=0, ndlenxy=0, dd0max=0;
941  dlen=d->DecayLength();
942  cosp=d->CosPointingAngle();
943  sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
944  maxdca=-9999.;
945  for(Int_t idau=0;idau<3;idau++) if(d->GetDCA(idau)>maxdca) maxdca=d->GetDCA(idau);
946  sigvert=d->GetSigmaVert();
947  ptmax=0;
948  dlenxy = d->DecayLengthXY();
949  ndlenxy=d->NormalizedDecayLengthXY();
950  cospxy=d->CosPointingAngleXY();
951  for(Int_t i=0; i<3; i++) {
952  if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
953  if(d->PtProng(i)<minPtDau) minPtDau=d->PtProng(i);
954  Double_t diffIP, errdiffIP;
955  d->Getd0MeasMinusExpProng(i,aod->GetMagneticField(),diffIP,errdiffIP);
956  Double_t normdd0= diffIP/errdiffIP;
957  if(i==0) dd0max=normdd0;
958  else if(TMath::Abs(normdd0)>TMath::Abs(dd0max)) dd0max=normdd0;
959  }
960  }
961  Double_t impparXY=d->ImpParXY()*10000.;
962  Double_t resSel=0;
963  if(fRDCutsAnalysis->GetIsSelectedCuts()) resSel=1;
964  if(passTopolAndPIDCuts) resSel=2;
965 
966  //for all THnSParse except for FD
967  Double_t arrayForSparse[kVarForSparse]={invMass,ptCand,impparXY,resSel,minPtDau,sigvert,cosp,cospxy,dlen,dlenxy,ndlenxy,dd0max};
968  //for THnSparse for FD
969  Double_t arrayForSparseFD[kVarForSparseFD]={invMass,ptCand,impparXY,resSel,minPtDau,sigvert,cosp,cospxy,dlen,dlenxy,ndlenxy,dd0max,ptB};
970  Double_t arrayForSparseTrue[kVarForSparseFD]={invMass,ptCand,trueImpParXY,resSel,minPtDau,sigvert,cosp,cospxy,dlen,dlenxy,ndlenxy,dd0max,ptB};
971  Double_t flagOrigin = 0;
972 
973  //Fill histos
974  index=GetHistoIndex(iPtBin);
975  fHistNEvents->Fill(9);
976  nSelectednopid++;
977  fPtVsMassNoPid->Fill(invMass,ptCand);
978  fMassHistNoPid[index]->Fill(invMass);
979  if(fDoImpPar){
980  fHistMassPtImpPar[0]->Fill(arrayForSparse);
981  }
982  if(passTopolAndPIDCuts){
983  fHistNEvents->Fill(10);
984  nSelected++;
985  fPtVsMass->Fill(invMass,ptCand);
986  fMassHist[index]->Fill(invMass);
987  if(d->GetCharge()>0) fMassHistPlus[index]->Fill(invMass);
988  else if(d->GetCharge()<0) fMassHistMinus[index]->Fill(invMass);
989  fPhiEtaCand->Fill(etaD,phiD);
990  if(TMath::Abs(invMass-1.8696)<0.05) fPhiEtaCandSigReg->Fill(etaD,phiD);
991  if(fCutsDistr){
992  fCosPHist[index]->Fill(cosp);
993  fDLenHist[index]->Fill(dlen);
994  fSumd02Hist[index]->Fill(sumD02);
995  fSigVertHist[index]->Fill(sigvert);
996  fPtMaxHist[index]->Fill(ptmax);
997  fPtKHist[index]->Fill(d->PtProng(1));
998  fPtpi1Hist[index]->Fill(d->PtProng(0));
999  fPtpi2Hist[index]->Fill(d->PtProng(2));
1000  fDCAHist[index]->Fill(maxdca);
1001  fDLxy[index]->Fill(ndlenxy);
1002  fCosxy[index]->Fill(cospxy);
1003  fCorreld0Kd0pi[0]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),d->Getd0Prong(2)*d->Getd0Prong(1));
1004  }
1005  }
1006 
1007  // fill ntuple
1008  if(fFillNtuple==1){
1009 
1010  Float_t tmp[31];
1011 
1012  tmp[0]=pdgCode;
1013  if(isFeeddown) tmp[0]+=5000.;
1014  tmp[1]=d->Px();
1015  tmp[2]=d->Py();
1016  tmp[3]=d->Pz();
1017  tmp[4]=d->Pt();
1018  tmp[5]=d->GetCharge();
1019  // tmp[5]=fRDCutsAnalysis->GetPIDBitMask(d);
1020  tmp[6]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(0));
1021  tmp[7]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(1));
1022  tmp[8]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(2));
1023  tmp[9]=d->PtProng(0);
1024  tmp[10]=d->PtProng(1);
1025  tmp[11]=d->PtProng(2);
1026  tmp[12]=d->PProng(0);
1027  tmp[13]=d->PProng(1);
1028  tmp[14]=d->PProng(2);
1029  tmp[15]=cosp;
1030  tmp[16]=cospxy;
1031  tmp[17]=dlen;
1032  tmp[18]=d->NormalizedDecayLength();
1033  tmp[19]=dlenxy;
1034  tmp[20]=ndlenxy;
1035  tmp[21]=d->InvMassDplus();
1036  tmp[22]=sigvert;
1037  tmp[23]=d->Getd0Prong(0);
1038  tmp[24]=d->Getd0Prong(1);
1039  tmp[25]=d->Getd0Prong(2);
1040  tmp[26]=maxdca;
1041  tmp[27]=ntracks;
1042  tmp[28]=fRDCutsAnalysis->GetCentrality(aod);
1043  tmp[29]=runNumber;
1044  tmp[30]=d->HasBadDaughters();
1045  fNtupleDplus->Fill(tmp);
1046  PostData(4,fNtupleDplus);
1047  }
1048 
1049  if(fFillNtuple==2 && passTopolAndPIDCuts){
1050  Float_t tmp[5];
1051  tmp[0]=pdgCode;
1052  if(isFeeddown) tmp[0]+=5000.;
1053  tmp[1]=d->Pt();
1054  tmp[2]=d->InvMassDplus();
1055  tmp[3]=impparXY;
1056  if(!fReadMC){flagOrigin=0;};
1057  if(fReadMC){
1058  if(isPrimary&&labDp>=0)flagOrigin=1;
1059  if(isFeeddown&&labDp>=0)flagOrigin=2;
1060  if(!(labDp>=0))flagOrigin=3;
1061 
1062  }
1063 
1064  tmp[4]=flagOrigin;
1065  fNtupleDplus->Fill(tmp);
1066  PostData(4,fNtupleDplus);
1067  }
1068 
1069 
1070  if(fReadMC){
1071  if(labDp>=0) {
1072  index=GetSignalHistoIndex(iPtBin);
1073  if(fDoImpPar){
1074  if(isPrimary) fHistMassPtImpPar[1]->Fill(arrayForSparse);
1075  else if(isFeeddown){
1076  fHistMassPtImpPar[2]->Fill(arrayForSparseFD);
1077  fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
1078  }
1079  }
1080  }else{
1081  index=GetBackgroundHistoIndex(iPtBin);
1082  if(fDoImpPar)fHistMassPtImpPar[4]->Fill(arrayForSparse);
1083  }
1084  fMassHistNoPid[index]->Fill(invMass);
1085  if(passTopolAndPIDCuts){
1086  if(fCutsDistr){
1087  Float_t fact=1.;
1088  Float_t factor[3]={1.,1.,1.};
1089  if(fUseStrangeness) fact=GetStrangenessWeights(d,arrayMC,factor);
1090  fCosPHist[index]->Fill(cosp,fact);
1091  fDLenHist[index]->Fill(dlen,fact);
1092  fDLxy[index]->Fill(ndlenxy);
1093  fCosxy[index]->Fill(cospxy);
1094  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];
1095  fSumd02Hist[index]->Fill(sumd02s);
1096  fSigVertHist[index]->Fill(sigvert,fact);
1097  fPtMaxHist[index]->Fill(ptmax,fact);
1098  fPtKHist[index]->Fill(d->PtProng(1),fact);
1099  fPtpi1Hist[index]->Fill(d->PtProng(0),fact);
1100  fPtpi2Hist[index]->Fill(d->PtProng(2),fact);
1101  fDCAHist[index]->Fill(maxdca,fact);
1102  fCorreld0Kd0pi[1]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),d->Getd0Prong(2)*d->Getd0Prong(1));
1103  }
1104  fMassHist[index]->Fill(invMass);
1105  if(d->GetCharge()>0) fMassHistPlus[index]->Fill(invMass);
1106  else if(d->GetCharge()<0) fMassHistMinus[index]->Fill(invMass);
1107  }
1108  }
1109  }
1110 
1111  if(recVtx)fRDCutsAnalysis->CleanOwnPrimaryVtx(d,aod,origownvtx);
1112 
1113  if(unsetvtx) d->UnsetOwnPrimaryVtx();
1114  }
1115  fCounter->StoreCandidates(aod,nSelectednopid,kTRUE);
1116  fCounter->StoreCandidates(aod,nSelected,kFALSE);
1117  }
1118  //start LS analysis
1119  if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
1120 
1121  PostData(1,fOutput);
1122  PostData(2,fListCuts);
1123  PostData(3,fCounter);
1124  return;
1125 }
1126 
1127 //________________________________________________________________________
1130 
1131  TString hisname;
1132  Int_t indexLS=0;
1133  Int_t index=0;
1134  Int_t nbins=GetNBinsHistos();
1135  for(Int_t i=0;i<fNPtBins;i++){
1136 
1137  index=GetHistoIndex(i);
1138  indexLS=GetLSHistoIndex(i);
1139 
1140  hisname.Form("hLSPt%d",i);
1141  fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1142  fMassHistLS[indexLS]->Sumw2();
1143 
1144  hisname.Form("hCosPAllPt%dLS",i);
1145  fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1146  fCosPHistLS[index]->Sumw2();
1147  hisname.Form("hDLenAllPt%dLS",i);
1148  fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1149  fDLenHistLS[index]->Sumw2();
1150  hisname.Form("hSumd02AllPt%dLS",i);
1151  fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1152  fSumd02HistLS[index]->Sumw2();
1153  hisname.Form("hSigVertAllPt%dLS",i);
1154  fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1155  fSigVertHistLS[index]->Sumw2();
1156  hisname.Form("hPtMaxAllPt%dLS",i);
1157  fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1158  fPtMaxHistLS[index]->Sumw2();
1159  hisname.Form("hDCAAllPt%dLS",i);
1160  fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1161  fDCAHistLS[index]->Sumw2();
1162 
1163  index=GetSignalHistoIndex(i);
1164  indexLS++;
1165 
1166  hisname.Form("hLSPt%dnw",i);
1167  fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1168  fMassHistLS[indexLS]->Sumw2();
1169 
1170  hisname.Form("hCosPSigPt%dLS",i);
1171  fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1172  fCosPHistLS[index]->Sumw2();
1173  hisname.Form("hDLenSigPt%dLS",i);
1174  fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1175  fDLenHistLS[index]->Sumw2();
1176  hisname.Form("hSumd02SigPt%dLS",i);
1177  fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1178  fSumd02HistLS[index]->Sumw2();
1179  hisname.Form("hSigVertSigPt%dLS",i);
1180  fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1181  fSigVertHistLS[index]->Sumw2();
1182  hisname.Form("hPtMaxSigPt%dLS",i);
1183  fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1184  fPtMaxHistLS[index]->Sumw2();
1185  hisname.Form("hDCASigPt%dLS",i);
1186  fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1187  fDCAHistLS[index]->Sumw2();
1188 
1189  index=GetBackgroundHistoIndex(i);
1190  indexLS++;
1191 
1192  hisname.Form("hLSPt%dLCntrip",i);
1193  fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1194  fMassHistLS[indexLS]->Sumw2();
1195 
1196  hisname.Form("hCosPBkgPt%dLS",i);
1197  fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1198  fCosPHistLS[index]->Sumw2();
1199  hisname.Form("hDLenBkgPt%dLS",i);
1200  fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1201  fDLenHistLS[index]->Sumw2();
1202  hisname.Form("hSumd02BkgPt%dLS",i);
1203  fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1204  fSumd02HistLS[index]->Sumw2();
1205  hisname.Form("hSigVertBkgPt%dLS",i);
1206  fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1207  fSigVertHistLS[index]->Sumw2();
1208  hisname.Form("hPtMaxBkgPt%dLS",i);
1209  fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1210  fPtMaxHistLS[index]->Sumw2();
1211  hisname.Form("hDCABkgPt%dLS",i);
1212  fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1213  fDCAHistLS[index]->Sumw2();
1214 
1215  indexLS++;
1216  hisname.Form("hLSPt%dLCntripsinglecut",i);
1217  fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1218  fMassHistLS[indexLS]->Sumw2();
1219 
1220  indexLS++;
1221  hisname.Form("hLSPt%dspc",i);
1222  fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1223  fMassHistLS[indexLS]->Sumw2();
1224  }
1225 
1226  for(Int_t i=0; i<3*fNPtBins; i++){
1227  fOutput->Add(fCosPHistLS[i]);
1228  fOutput->Add(fDLenHistLS[i]);
1229  fOutput->Add(fSumd02HistLS[i]);
1230  fOutput->Add(fSigVertHistLS[i]);
1231  fOutput->Add(fPtMaxHistLS[i]);
1232  fOutput->Add(fDCAHistLS[i]);
1233  }
1234  for(Int_t i=0; i<5*fNPtBins; i++){
1235  fOutput->Add(fMassHistLS[i]);
1236  }
1237 }
1238 
1239 //________________________________________________________________________
1242 
1243  Int_t nmassbins=GetNBinsHistos();
1244  Double_t maxmult;
1245  if(fSystem==1) maxmult=5000;
1246  else maxmult=200;
1247 
1248  Int_t nptbins=80;
1249  Double_t ptmin=0.;
1250  Double_t ptmax=40.;
1251 
1252  Int_t nptmindaubins=10;
1253  Double_t minptmindau=0.2;
1254  Double_t maxptmindau=1.2;
1255 
1256  Int_t nsigvertbins=25;
1257  Double_t minsigvert=0.010;
1258  Double_t maxsigvert=0.035;
1259 
1260  Int_t nselbins=2;
1261  Double_t minsel=0.5;
1262  Double_t maxsel=2.5;
1263 
1264  Int_t ncospbins=100;
1265  Double_t mincosp=0.90;
1266  Double_t maxcosp=1.;
1267 
1268  Int_t ncospxybins=30;
1269  Double_t mincospxy=0.97;
1270  Double_t maxcospxy=1.;
1271 
1272  Int_t ndeclbins=70;
1273  Double_t mindecl=0.;
1274  Double_t maxdecl=0.7;
1275 
1276  Int_t ndeclxybins=70;
1277  Double_t mindeclxy=0.;
1278  Double_t maxdeclxy=0.7;
1279 
1280  Int_t nnormdlbins=30;
1281  Double_t minnormdl=0.;
1282  Double_t maxnormdl=30.;
1283 
1284  Int_t nd0d0expbins=40;
1285  Double_t mind0d0=-10.;
1286  Double_t maxd0d0=10.;
1287 
1288 
1289  //dimensions for THnSparse which are NOT for BFeed
1290  TString axTit[kVarForSparse]={"M_{K#pi#pi} (GeV/c^{2})","p_{T} (GeV/c)","Imp Par (#mum)","passTopolPID","min. daughter p_{T} (GeV/c)",
1291  "sigmaVertex","cos(#theta_{P})","cos(#theta_{P}^{xy})","decL (cm)","decL XY (cm)","Norm decL XY","Norm max d0-d0exp"};
1292 
1293  Int_t nbins[kVarForSparse]={nmassbins,nptbins,fNImpParBins,nselbins,nptmindaubins,nsigvertbins,ncospbins,ncospxybins,ndeclbins,ndeclxybins,nnormdlbins,nd0d0expbins};
1294  Double_t xmin[kVarForSparse]={fLowmasslimit,ptmin,fLowerImpPar,minsel,minptmindau,minsigvert,mincosp,mincospxy,mindecl,mindeclxy,minnormdl,mind0d0};
1295  Double_t xmax[kVarForSparse]={fUpmasslimit,ptmax,fHigherImpPar,maxsel,maxptmindau,maxsigvert,maxcosp,maxcospxy,maxdecl,maxdeclxy,maxnormdl,maxd0d0};
1296 
1297  //dimensions for THnSparse for BFeed
1298  Int_t nbinsFD[kVarForSparseFD]={nmassbins,nptbins,fNImpParBins,nselbins,nptmindaubins,nsigvertbins,ncospbins,ncospxybins,ndeclbins,ndeclbins,nnormdlbins,nd0d0expbins,84};
1299  Double_t xminFD[kVarForSparseFD]={fLowmasslimit,ptmin,fLowerImpPar,minsel,minptmindau,minsigvert,mincosp,mincospxy,mindecl,mindeclxy,minnormdl,mind0d0,-2};
1300  Double_t xmaxFD[kVarForSparseFD]={fUpmasslimit,ptmax,fHigherImpPar,maxsel,maxptmindau,maxsigvert,maxcosp,maxcospxy,maxdecl,maxdeclxy,maxnormdl,maxd0d0,40};
1301 
1302  //mass, pt, imppar, cosPoinXY, decLXY, norm decLXY (for BFeed also ptB)
1303  //mass, pt, imppar, cosPoinXY, decLXY, norm decLXY (for BFeed also ptB)
1304  fHistMassPtImpPar[0]=new THnSparseF("hMassPtImpParAll",
1305  "Mass vs. pt vs.imppar - All",
1306  kVarForSparse,nbins,xmin,xmax);
1307  fHistMassPtImpPar[1]=new THnSparseF("hMassPtImpParPrompt",
1308  "Mass vs. pt vs.imppar - promptD",
1309  kVarForSparse,nbins,xmin,xmax);
1310  fHistMassPtImpPar[2]=new THnSparseF("hMassPtImpParBfeed",
1311  "Mass vs. pt vs.imppar - DfromB",
1312  kVarForSparseFD,nbinsFD,xminFD,xmaxFD);
1313  fHistMassPtImpPar[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1314  "Mass vs. pt vs.true imppar -DfromB",
1315  kVarForSparseFD,nbinsFD,xminFD,xmaxFD);
1316  fHistMassPtImpPar[4]=new THnSparseF("hMassPtImpParBkg",
1317  "Mass vs. pt vs.imppar - backgr.",
1318  kVarForSparse,nbins,xmin,xmax);
1319 
1320 
1321  for(Int_t i=0; i<5; i++){
1322  for(Int_t iax=0; iax<kVarForSparse; iax++) fHistMassPtImpPar[i]->GetAxis(iax)->SetTitle(axTit[iax].Data());
1323  if(i == 2 || i == 3) fHistMassPtImpPar[i]->GetAxis(kVarForSparseFD-1)->SetTitle("p_{T}^{B} (GeV/c)");
1324  fOutput->Add(fHistMassPtImpPar[i]);
1325  }
1326 }
1327 
1328 //________________________________________________________________________
1331 
1332  const Int_t nVarPrompt = 2;
1333  const Int_t nVarFD = 3;
1334 
1335  Int_t nbinsPrompt[nVarPrompt]={200,100};
1336  Int_t nbinsFD[nVarFD]={200,100,200};
1337 
1338  Double_t xminPrompt[nVarPrompt] = {0.,-1.};
1339  Double_t xmaxPrompt[nVarPrompt] = {40.,1.};
1340 
1341  Double_t xminFD[nVarFD] = {0.,-1.,0.};
1342  Double_t xmaxFD[nVarFD] = {40.,1.,40.};
1343 
1344  //pt, y
1345  fMCAccPrompt = new THnSparseF("hMCAccPrompt","kStepMCAcceptance pt vs. y - promptD",nVarPrompt,nbinsPrompt,xminPrompt,xmaxPrompt);
1346  fMCAccPrompt->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
1347  fMCAccPrompt->GetAxis(1)->SetTitle("y");
1348 
1349  //pt,y,ptB
1350  fMCAccBFeed = new THnSparseF("hMCAccBFeed","kStepMCAcceptance pt vs. y vs. ptB - DfromB",nVarFD,nbinsFD,xminFD,xmaxFD);
1351  fMCAccBFeed->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
1352  fMCAccBFeed->GetAxis(1)->SetTitle("y");
1353  fMCAccBFeed->GetAxis(2)->SetTitle("p_{T}^{B} (GeV/c)");
1354 
1355  fOutput->Add(fMCAccPrompt);
1356  fOutput->Add(fMCAccBFeed);
1357 }
1358 
1359 //________________________________________________________________________
1360 void AliAnalysisTaskSEDplus::FillMCAcceptanceHistos(TClonesArray *arrayMC, AliAODMCHeader *mcHeader){
1362 
1363  const Int_t nProng = 3;
1364 
1365  Double_t zMCVertex = mcHeader->GetVtxZ(); //vertex MC
1366 
1367  for(Int_t iPart=0; iPart<arrayMC->GetEntriesFast(); iPart++){
1368  AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(arrayMC->At(iPart));
1369  if (TMath::Abs(mcPart->GetPdgCode()) == 411){
1370 
1371  Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,mcPart,fUseQuarkTagInKine);//Prompt = 4, FeedDown = 5
1372 
1373  Int_t deca = 0;
1374  Bool_t isGoodDecay=kFALSE;
1375  Int_t labDau[4]={-1,-1,-1,-1};
1376  Bool_t isInAcc = kFALSE;
1377  Bool_t isFidAcc = kFALSE;
1378 
1379  deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,mcPart,labDau);
1380  if(deca > 0) isGoodDecay=kTRUE;
1381 
1382  if(labDau[0]==-1){
1383  continue; //protection against unfilled array of labels
1384  }
1385 
1386  isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(mcPart->Pt(),mcPart->Y());
1387  isInAcc=CheckAcc(arrayMC,nProng,labDau);
1388 
1389  if(isGoodDecay && TMath::Abs(zMCVertex) < fRDCutsAnalysis->GetMaxVtxZ() && isFidAcc && isInAcc) {
1390  //for prompt
1391  if(orig == 4){
1392  //fill histo for prompt
1393  Double_t arrayMCprompt[2] = {mcPart->Pt(),mcPart->Y()};
1394  fMCAccPrompt->Fill(arrayMCprompt);
1395  }
1396  //for FD
1397  else if(orig == 5){
1398  Double_t ptB = AliVertexingHFUtils::GetBeautyMotherPt(arrayMC,mcPart);
1399  //fill histo for FD
1400  Double_t arrayMCFD[3] = {mcPart->Pt(),mcPart->Y(),ptB};
1401  fMCAccBFeed->Fill(arrayMCFD);
1402  }
1403  else
1404  continue;
1405  }
1406  }
1407  }
1408 }
1409 
1410 //________________________________________________________________________
1411 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
1412 {
1414  //
1415  if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
1416 
1417  fOutput = dynamic_cast<TList*> (GetOutputData(1));
1418  if (!fOutput) {
1419  printf("ERROR: fOutput not available\n");
1420  return;
1421  }
1422 
1423  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
1424  if(fHistNEvents){
1425  printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
1426  }else{
1427  printf("ERROR: fHistNEvents not available\n");
1428  return;
1429  }
1430 
1431  return;
1432 }
1433 //_________________________________________________________________________________________________
1434 Float_t AliAnalysisTaskSEDplus::GetTrueImpactParameter(const AliAODMCHeader *mcHeader, TClonesArray* arrayMC, const AliAODMCParticle *partDp) const {
1436 
1437  Double_t vtxTrue[3];
1438  mcHeader->GetVertex(vtxTrue);
1439  Double_t origD[3];
1440  partDp->XvYvZv(origD);
1441  Short_t charge=partDp->Charge();
1442  Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
1443  for(Int_t iDau=0; iDau<3; iDau++){
1444  pXdauTrue[iDau]=0.;
1445  pYdauTrue[iDau]=0.;
1446  pZdauTrue[iDau]=0.;
1447  }
1448 
1449  Int_t nDau=partDp->GetNDaughters();
1450  Int_t labelFirstDau = partDp->GetDaughter(0);
1451  if(nDau==3){
1452  for(Int_t iDau=0; iDau<3; iDau++){
1453  Int_t ind = labelFirstDau+iDau;
1454  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1455  if(!part){
1456  AliError("Daughter particle not found in MC array");
1457  return 99999.;
1458  }
1459  pXdauTrue[iDau]=part->Px();
1460  pYdauTrue[iDau]=part->Py();
1461  pZdauTrue[iDau]=part->Pz();
1462  }
1463  }else if(nDau==2){
1464  Int_t theDau=0;
1465  for(Int_t iDau=0; iDau<2; iDau++){
1466  Int_t ind = labelFirstDau+iDau;
1467  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1468  if(!part){
1469  AliError("Daughter particle not found in MC array");
1470  return 99999.;
1471  }
1472  Int_t pdgCode=TMath::Abs(part->GetPdgCode());
1473  if(pdgCode==211 || pdgCode==321){
1474  pXdauTrue[theDau]=part->Px();
1475  pYdauTrue[theDau]=part->Py();
1476  pZdauTrue[theDau]=part->Pz();
1477  ++theDau;
1478  }else{
1479  Int_t nDauRes=part->GetNDaughters();
1480  if(nDauRes==2){
1481  Int_t labelFirstDauRes = part->GetDaughter(0);
1482  for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
1483  Int_t indDR = labelFirstDauRes+iDauRes;
1484  AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
1485  if(!partDR){
1486  AliError("Daughter particle not found in MC array");
1487  return 99999.;
1488  }
1489 
1490  Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
1491  if(pdgCodeDR==211 || pdgCodeDR==321){
1492  pXdauTrue[theDau]=partDR->Px();
1493  pYdauTrue[theDau]=partDR->Py();
1494  pZdauTrue[theDau]=partDR->Pz();
1495  ++theDau;
1496  }
1497  }
1498  }
1499  }
1500  }
1501  if(theDau!=3){
1502  AliError("Wrong number of decay prongs");
1503  return 99999.;
1504  }
1505  }
1506 
1507  Double_t d0dummy[3]={0.,0.,0.};
1508  AliAODRecoDecayHF aodDplusMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
1509  return aodDplusMC.ImpParXY();
1510 
1511 }
1512 //_________________________________________________________________________________________________
1513 Float_t AliAnalysisTaskSEDplus::GetStrangenessWeights(const AliAODRecoDecayHF3Prong* d, TClonesArray* arrayMC, Float_t factor[3]) const {
1515 
1516  for(Int_t iprong=0;iprong<3;iprong++){
1517  factor[iprong]=1;
1518  AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
1519  Int_t labd= trad->GetLabel();
1520  if(labd>=0){
1521  AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
1522  if(dau){
1523  Int_t labm = dau->GetMother();
1524  if(labm>=0){
1525  AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
1526  if(mot){
1527  if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
1528  if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
1529  else factor[iprong]=1./.6;
1530  // fNentries->Fill(12);
1531  }
1532  if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
1533  factor[iprong]=1./0.25;
1534  // fNentries->Fill(13);
1535  }//if 3122
1536  }//if(mot)
1537  }//if labm>0
1538  }//if(dau)
1539  }//if labd>=0
1540  }//prong loop
1541 
1542  Float_t fact=1.;
1543  for(Int_t k=0;k<3;k++)fact=fact*factor[k];
1544  return fact;
1545 
1546 }
1547 
1548 //_________________________________________________________________
1549 Bool_t AliAnalysisTaskSEDplus::CheckAcc(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
1551  for (Int_t iProng = 0; iProng<nProng; iProng++){
1552  AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
1553  if(!mcPartDaughter) return kFALSE;
1554  Double_t eta = mcPartDaughter->Eta();
1555  Double_t pt = mcPartDaughter->Pt();
1556  if (TMath::Abs(eta) > 0.9 || pt < 0.1) return kFALSE;
1557  }
1558  return kTRUE;
1559 }
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:333
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:293
Double_t CosPointingAngleXY() const
TH1F * fMassHistNoPid[3 *kMaxPtBins]
!hist. for inv mass (w/o PID)
TList * fListCuts
width of one bin in output histos
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:261
TH1F * fSumd02HistLS[3 *kMaxPtBins]
!hist. for LS cuts variable 3 (topol+PID)
Bool_t IsEventSelected(AliVEvent *event)
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
void StoreEvent(AliVEvent *, AliRDHFCuts *, Bool_t mc=kFALSE, Int_t multiplicity=-9999)
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:263
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)