AliPhysics  v5-07-15-01 (b3d7633)
 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 ",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  Int_t passTopolAndPIDCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
873 
874  if(!(vHF->FillRecoCand(aod,d))) { //Fill the data members of the candidate only if they are empty.
875  fHistNEvents->Fill(12); //monitor how often this fails
876  continue;
877  }
878 
879  Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
880 
881  if(!fRDCutsAnalysis->GetIsSelectedCuts()) continue;
882 
883  Double_t etaD=d->Eta();
884  Double_t phiD=d->Phi();
885  if(fEtaSelection!=0){
886  if(fEtaSelection==1 && etaD<0) continue;
887  if(fEtaSelection==-1 && etaD>0) continue;
888  }
889 
890  Bool_t unsetvtx=kFALSE;
891  if(!d->GetOwnPrimaryVtx()){
892  d->SetOwnPrimaryVtx(vtx1);
893  unsetvtx=kTRUE;
894  }
895 
896  Double_t ptCand = d->Pt();
897  Int_t iPtBin = fRDCutsAnalysis->PtBin(ptCand);
898 
899  Bool_t recVtx=kFALSE;
900  AliAODVertex *origownvtx=0x0;
902  if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
903  if(fRDCutsAnalysis->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
904  else fRDCutsAnalysis->CleanOwnPrimaryVtx(d,aod,origownvtx);
905  }
906 
907  Int_t labDp=-1;
908  Bool_t isPrimary=kFALSE;
909  Bool_t isFeeddown=kFALSE;
910  Float_t pdgCode=-2;
911  Float_t trueImpParXY=0.;
912  Double_t ptB=-1.5;
913  if(fReadMC){
914  labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
915  if(labDp>=0){
916  AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
917  Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,partDp,fUseQuarkTagInKine);//Prompt = 4, FeedDown = 5
918  pdgCode=TMath::Abs(partDp->GetPdgCode());
919  if(orig==4){
920  isPrimary=kTRUE;
921  isFeeddown=kFALSE;
922  }else if(orig==5){
923  isPrimary=kFALSE;
924  isFeeddown=kTRUE;
925  trueImpParXY=GetTrueImpactParameter(mcHeader,arrayMC,partDp)*10000.;
926  ptB=AliVertexingHFUtils::GetBeautyMotherPt(arrayMC,partDp);
927  }else{
928  pdgCode=-3;
929  }
930  }else{
931  pdgCode=-1;
932  }
933  }
934 
935  Double_t invMass=d->InvMassDplus();
936  Double_t rapid=d->YDplus();
937  fYVsPtNoPid->Fill(ptCand,rapid,invMass);
938  if(passTopolAndPIDCuts) {
939  fYVsPt->Fill(ptCand,rapid,invMass);
940  nOS++;
941  }
942  if(fReadMC && labDp>=0){
943  fYVsPtSigNoPid->Fill(ptCand,rapid, invMass);
944  if(passTopolAndPIDCuts)fYVsPtSig->Fill(ptCand,rapid, invMass);
945  }
946 
947  Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
948  if(isFidAcc){
949 
950  Double_t minPtDau=999.,ptmax=0,maxdca=0,sigvert=0,sumD02=0;
951  Double_t dlen=0,cosp=0,dlenxy=0,cospxy=0, ndlenxy=0, dd0max=0;
953  dlen=d->DecayLength();
954  cosp=d->CosPointingAngle();
955  sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
956  maxdca=-9999.;
957  for(Int_t idau=0;idau<3;idau++) if(d->GetDCA(idau)>maxdca) maxdca=d->GetDCA(idau);
958  sigvert=d->GetSigmaVert();
959  ptmax=0;
960  dlenxy = d->DecayLengthXY();
961  ndlenxy=d->NormalizedDecayLengthXY();
962  cospxy=d->CosPointingAngleXY();
963  for(Int_t i=0; i<3; i++) {
964  if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
965  if(d->PtProng(i)<minPtDau) minPtDau=d->PtProng(i);
966  Double_t diffIP, errdiffIP;
967  d->Getd0MeasMinusExpProng(i,aod->GetMagneticField(),diffIP,errdiffIP);
968  Double_t normdd0= diffIP/errdiffIP;
969  if(i==0) dd0max=normdd0;
970  else if(TMath::Abs(normdd0)>TMath::Abs(dd0max)) dd0max=normdd0;
971  }
972  }
973  Double_t impparXY=d->ImpParXY()*10000.;
974  Double_t resSel=0;
975  if(fRDCutsAnalysis->GetIsSelectedCuts()) resSel=1;
976  if(passTopolAndPIDCuts) resSel=2;
977 
978  //for all THnSParse except for FD
979  Double_t arrayForSparse[kVarForSparse]={invMass,ptCand,impparXY,resSel,minPtDau,sigvert,cosp,cospxy,dlen,dlenxy,ndlenxy,dd0max};
980  //for THnSparse for FD
981  Double_t arrayForSparseFD[kVarForSparseFD]={invMass,ptCand,impparXY,resSel,minPtDau,sigvert,cosp,cospxy,dlen,dlenxy,ndlenxy,dd0max,ptB};
982  Double_t arrayForSparseTrue[kVarForSparseFD]={invMass,ptCand,trueImpParXY,resSel,minPtDau,sigvert,cosp,cospxy,dlen,dlenxy,ndlenxy,dd0max,ptB};
983  Double_t flagOrigin = 0;
984 
985  //Fill histos
986  index=GetHistoIndex(iPtBin);
987  fHistNEvents->Fill(9);
988  nSelectednopid++;
989  fPtVsMassNoPid->Fill(invMass,ptCand);
990  fMassHistNoPid[index]->Fill(invMass);
991  if(fDoImpPar){
992  fHistMassPtImpPar[0]->Fill(arrayForSparse);
993  }
994  if(passTopolAndPIDCuts){
995  fHistNEvents->Fill(10);
996  nSelected++;
997  fPtVsMass->Fill(invMass,ptCand);
998  fMassHist[index]->Fill(invMass);
999  if(d->GetCharge()>0) fMassHistPlus[index]->Fill(invMass);
1000  else if(d->GetCharge()<0) fMassHistMinus[index]->Fill(invMass);
1001  fPhiEtaCand->Fill(etaD,phiD);
1002  if(TMath::Abs(invMass-1.8696)<0.05) fPhiEtaCandSigReg->Fill(etaD,phiD);
1003  if(fCutsDistr){
1004  fCosPHist[index]->Fill(cosp);
1005  fDLenHist[index]->Fill(dlen);
1006  fSumd02Hist[index]->Fill(sumD02);
1007  fSigVertHist[index]->Fill(sigvert);
1008  fPtMaxHist[index]->Fill(ptmax);
1009  fPtKHist[index]->Fill(d->PtProng(1));
1010  fPtpi1Hist[index]->Fill(d->PtProng(0));
1011  fPtpi2Hist[index]->Fill(d->PtProng(2));
1012  fDCAHist[index]->Fill(maxdca);
1013  fDLxy[index]->Fill(ndlenxy);
1014  fCosxy[index]->Fill(cospxy);
1015  fCorreld0Kd0pi[0]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),d->Getd0Prong(2)*d->Getd0Prong(1));
1016  }
1017  }
1018 
1019  // fill ntuple
1020  if(fFillNtuple==1){
1021 
1022  Float_t tmp[31];
1023 
1024  tmp[0]=pdgCode;
1025  if(isFeeddown) tmp[0]+=5000.;
1026  tmp[1]=d->Px();
1027  tmp[2]=d->Py();
1028  tmp[3]=d->Pz();
1029  tmp[4]=d->Pt();
1030  tmp[5]=d->GetCharge();
1031  // tmp[5]=fRDCutsAnalysis->GetPIDBitMask(d);
1032  tmp[6]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(0));
1033  tmp[7]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(1));
1034  tmp[8]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(2));
1035  tmp[9]=d->PtProng(0);
1036  tmp[10]=d->PtProng(1);
1037  tmp[11]=d->PtProng(2);
1038  tmp[12]=d->PProng(0);
1039  tmp[13]=d->PProng(1);
1040  tmp[14]=d->PProng(2);
1041  tmp[15]=cosp;
1042  tmp[16]=cospxy;
1043  tmp[17]=dlen;
1044  tmp[18]=d->NormalizedDecayLength();
1045  tmp[19]=dlenxy;
1046  tmp[20]=ndlenxy;
1047  tmp[21]=d->InvMassDplus();
1048  tmp[22]=sigvert;
1049  tmp[23]=d->Getd0Prong(0);
1050  tmp[24]=d->Getd0Prong(1);
1051  tmp[25]=d->Getd0Prong(2);
1052  tmp[26]=maxdca;
1053  tmp[27]=ntracks;
1054  tmp[28]=fRDCutsAnalysis->GetCentrality(aod);
1055  tmp[29]=runNumber;
1056  tmp[30]=d->HasBadDaughters();
1057  fNtupleDplus->Fill(tmp);
1058  PostData(4,fNtupleDplus);
1059  }
1060 
1061  if(fFillNtuple==2 && passTopolAndPIDCuts){
1062  Float_t tmp[5];
1063  tmp[0]=pdgCode;
1064  if(isFeeddown) tmp[0]+=5000.;
1065  tmp[1]=d->Pt();
1066  tmp[2]=d->InvMassDplus();
1067  tmp[3]=impparXY;
1068  if(!fReadMC){flagOrigin=0;};
1069  if(fReadMC){
1070  if(isPrimary&&labDp>=0)flagOrigin=1;
1071  if(isFeeddown&&labDp>=0)flagOrigin=2;
1072  if(!(labDp>=0))flagOrigin=3;
1073 
1074  }
1075 
1076  tmp[4]=flagOrigin;
1077  fNtupleDplus->Fill(tmp);
1078  PostData(4,fNtupleDplus);
1079  }
1080 
1081 
1082  if(fReadMC){
1083  if(labDp>=0) {
1084  index=GetSignalHistoIndex(iPtBin);
1085  if(fDoImpPar){
1086  if(isPrimary) fHistMassPtImpPar[1]->Fill(arrayForSparse);
1087  else if(isFeeddown){
1088  fHistMassPtImpPar[2]->Fill(arrayForSparseFD);
1089  fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
1090  }
1091  }
1092  }else{
1093  index=GetBackgroundHistoIndex(iPtBin);
1094  if(fDoImpPar)fHistMassPtImpPar[4]->Fill(arrayForSparse);
1095  }
1096  fMassHistNoPid[index]->Fill(invMass);
1097  if(passTopolAndPIDCuts){
1098  if(fCutsDistr){
1099  Float_t fact=1.;
1100  Float_t factor[3]={1.,1.,1.};
1101  if(fUseStrangeness) fact=GetStrangenessWeights(d,arrayMC,factor);
1102  fCosPHist[index]->Fill(cosp,fact);
1103  fDLenHist[index]->Fill(dlen,fact);
1104  fDLxy[index]->Fill(ndlenxy);
1105  fCosxy[index]->Fill(cospxy);
1106  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];
1107  fSumd02Hist[index]->Fill(sumd02s);
1108  fSigVertHist[index]->Fill(sigvert,fact);
1109  fPtMaxHist[index]->Fill(ptmax,fact);
1110  fPtKHist[index]->Fill(d->PtProng(1),fact);
1111  fPtpi1Hist[index]->Fill(d->PtProng(0),fact);
1112  fPtpi2Hist[index]->Fill(d->PtProng(2),fact);
1113  fDCAHist[index]->Fill(maxdca,fact);
1114  fCorreld0Kd0pi[1]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),d->Getd0Prong(2)*d->Getd0Prong(1));
1115  }
1116  fMassHist[index]->Fill(invMass);
1117  if(d->GetCharge()>0) fMassHistPlus[index]->Fill(invMass);
1118  else if(d->GetCharge()<0) fMassHistMinus[index]->Fill(invMass);
1119  }
1120  }
1121  }
1122 
1123  if(recVtx)fRDCutsAnalysis->CleanOwnPrimaryVtx(d,aod,origownvtx);
1124 
1125  if(unsetvtx) d->UnsetOwnPrimaryVtx();
1126  }
1127  fCounter->StoreCandidates(aod,nSelectednopid,kTRUE);
1128  fCounter->StoreCandidates(aod,nSelected,kFALSE);
1129  }
1130  delete vHF;
1131  //start LS analysis
1132  if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
1133 
1134  PostData(1,fOutput);
1135  PostData(2,fListCuts);
1136  PostData(3,fCounter);
1137  return;
1138 }
1139 
1140 //________________________________________________________________________
1143 
1144  TString hisname;
1145  Int_t indexLS=0;
1146  Int_t index=0;
1147  Int_t nbins=GetNBinsHistos();
1148  for(Int_t i=0;i<fNPtBins;i++){
1149 
1150  index=GetHistoIndex(i);
1151  indexLS=GetLSHistoIndex(i);
1152 
1153  hisname.Form("hLSPt%d",i);
1154  fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1155  fMassHistLS[indexLS]->Sumw2();
1156 
1157  hisname.Form("hCosPAllPt%dLS",i);
1158  fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1159  fCosPHistLS[index]->Sumw2();
1160  hisname.Form("hDLenAllPt%dLS",i);
1161  fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1162  fDLenHistLS[index]->Sumw2();
1163  hisname.Form("hSumd02AllPt%dLS",i);
1164  fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1165  fSumd02HistLS[index]->Sumw2();
1166  hisname.Form("hSigVertAllPt%dLS",i);
1167  fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1168  fSigVertHistLS[index]->Sumw2();
1169  hisname.Form("hPtMaxAllPt%dLS",i);
1170  fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1171  fPtMaxHistLS[index]->Sumw2();
1172  hisname.Form("hDCAAllPt%dLS",i);
1173  fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1174  fDCAHistLS[index]->Sumw2();
1175 
1176  index=GetSignalHistoIndex(i);
1177  indexLS++;
1178 
1179  hisname.Form("hLSPt%dnw",i);
1180  fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1181  fMassHistLS[indexLS]->Sumw2();
1182 
1183  hisname.Form("hCosPSigPt%dLS",i);
1184  fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1185  fCosPHistLS[index]->Sumw2();
1186  hisname.Form("hDLenSigPt%dLS",i);
1187  fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1188  fDLenHistLS[index]->Sumw2();
1189  hisname.Form("hSumd02SigPt%dLS",i);
1190  fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1191  fSumd02HistLS[index]->Sumw2();
1192  hisname.Form("hSigVertSigPt%dLS",i);
1193  fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1194  fSigVertHistLS[index]->Sumw2();
1195  hisname.Form("hPtMaxSigPt%dLS",i);
1196  fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1197  fPtMaxHistLS[index]->Sumw2();
1198  hisname.Form("hDCASigPt%dLS",i);
1199  fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1200  fDCAHistLS[index]->Sumw2();
1201 
1202  index=GetBackgroundHistoIndex(i);
1203  indexLS++;
1204 
1205  hisname.Form("hLSPt%dLCntrip",i);
1206  fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1207  fMassHistLS[indexLS]->Sumw2();
1208 
1209  hisname.Form("hCosPBkgPt%dLS",i);
1210  fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1211  fCosPHistLS[index]->Sumw2();
1212  hisname.Form("hDLenBkgPt%dLS",i);
1213  fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1214  fDLenHistLS[index]->Sumw2();
1215  hisname.Form("hSumd02BkgPt%dLS",i);
1216  fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1217  fSumd02HistLS[index]->Sumw2();
1218  hisname.Form("hSigVertBkgPt%dLS",i);
1219  fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1220  fSigVertHistLS[index]->Sumw2();
1221  hisname.Form("hPtMaxBkgPt%dLS",i);
1222  fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1223  fPtMaxHistLS[index]->Sumw2();
1224  hisname.Form("hDCABkgPt%dLS",i);
1225  fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1226  fDCAHistLS[index]->Sumw2();
1227 
1228  indexLS++;
1229  hisname.Form("hLSPt%dLCntripsinglecut",i);
1230  fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1231  fMassHistLS[indexLS]->Sumw2();
1232 
1233  indexLS++;
1234  hisname.Form("hLSPt%dspc",i);
1235  fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1236  fMassHistLS[indexLS]->Sumw2();
1237  }
1238 
1239  for(Int_t i=0; i<3*fNPtBins; i++){
1240  fOutput->Add(fCosPHistLS[i]);
1241  fOutput->Add(fDLenHistLS[i]);
1242  fOutput->Add(fSumd02HistLS[i]);
1243  fOutput->Add(fSigVertHistLS[i]);
1244  fOutput->Add(fPtMaxHistLS[i]);
1245  fOutput->Add(fDCAHistLS[i]);
1246  }
1247  for(Int_t i=0; i<5*fNPtBins; i++){
1248  fOutput->Add(fMassHistLS[i]);
1249  }
1250 }
1251 
1252 //________________________________________________________________________
1255 
1256  Int_t nmassbins=GetNBinsHistos();
1257  Double_t maxmult;
1258  if(fSystem==1) maxmult=5000;
1259  else maxmult=200;
1260 
1261  Int_t nptbins=80;
1262  Double_t ptmin=0.;
1263  Double_t ptmax=40.;
1264 
1265  Int_t nptmindaubins=10;
1266  Double_t minptmindau=0.2;
1267  Double_t maxptmindau=1.2;
1268 
1269  Int_t nsigvertbins=25;
1270  Double_t minsigvert=0.010;
1271  Double_t maxsigvert=0.035;
1272 
1273  Int_t nselbins=2;
1274  Double_t minsel=0.5;
1275  Double_t maxsel=2.5;
1276 
1277  Int_t ncospbins=100;
1278  Double_t mincosp=0.90;
1279  Double_t maxcosp=1.;
1280 
1281  Int_t ncospxybins=30;
1282  Double_t mincospxy=0.97;
1283  Double_t maxcospxy=1.;
1284 
1285  Int_t ndeclbins=70;
1286  Double_t mindecl=0.;
1287  Double_t maxdecl=0.7;
1288 
1289  Int_t ndeclxybins=70;
1290  Double_t mindeclxy=0.;
1291  Double_t maxdeclxy=0.7;
1292 
1293  Int_t nnormdlbins=30;
1294  Double_t minnormdl=0.;
1295  Double_t maxnormdl=30.;
1296 
1297  Int_t nd0d0expbins=40;
1298  Double_t mind0d0=-10.;
1299  Double_t maxd0d0=10.;
1300 
1301 
1302  //dimensions for THnSparse which are NOT for BFeed
1303  TString axTit[kVarForSparse]={"M_{K#pi#pi} (GeV/c^{2})","p_{T} (GeV/c)","Imp Par (#mum)","passTopolPID","min. daughter p_{T} (GeV/c)",
1304  "sigmaVertex","cos(#theta_{P})","cos(#theta_{P}^{xy})","decL (cm)","decL XY (cm)","Norm decL XY","Norm max d0-d0exp"};
1305 
1306  Int_t nbins[kVarForSparse]={nmassbins,nptbins,fNImpParBins,nselbins,nptmindaubins,nsigvertbins,ncospbins,ncospxybins,ndeclbins,ndeclxybins,nnormdlbins,nd0d0expbins};
1307  Double_t xmin[kVarForSparse]={fLowmasslimit,ptmin,fLowerImpPar,minsel,minptmindau,minsigvert,mincosp,mincospxy,mindecl,mindeclxy,minnormdl,mind0d0};
1308  Double_t xmax[kVarForSparse]={fUpmasslimit,ptmax,fHigherImpPar,maxsel,maxptmindau,maxsigvert,maxcosp,maxcospxy,maxdecl,maxdeclxy,maxnormdl,maxd0d0};
1309 
1310  //dimensions for THnSparse for BFeed
1311  Int_t nbinsFD[kVarForSparseFD]={nmassbins,nptbins,fNImpParBins,nselbins,nptmindaubins,nsigvertbins,ncospbins,ncospxybins,ndeclbins,ndeclbins,nnormdlbins,nd0d0expbins,84};
1312  Double_t xminFD[kVarForSparseFD]={fLowmasslimit,ptmin,fLowerImpPar,minsel,minptmindau,minsigvert,mincosp,mincospxy,mindecl,mindeclxy,minnormdl,mind0d0,-2};
1313  Double_t xmaxFD[kVarForSparseFD]={fUpmasslimit,ptmax,fHigherImpPar,maxsel,maxptmindau,maxsigvert,maxcosp,maxcospxy,maxdecl,maxdeclxy,maxnormdl,maxd0d0,40};
1314 
1315  //mass, pt, imppar, cosPoinXY, decLXY, norm decLXY (for BFeed also ptB)
1316  //mass, pt, imppar, cosPoinXY, decLXY, norm decLXY (for BFeed also ptB)
1317  fHistMassPtImpPar[0]=new THnSparseF("hMassPtImpParAll",
1318  "Mass vs. pt vs.imppar - All",
1319  kVarForSparse,nbins,xmin,xmax);
1320  fHistMassPtImpPar[1]=new THnSparseF("hMassPtImpParPrompt",
1321  "Mass vs. pt vs.imppar - promptD",
1322  kVarForSparse,nbins,xmin,xmax);
1323  fHistMassPtImpPar[2]=new THnSparseF("hMassPtImpParBfeed",
1324  "Mass vs. pt vs.imppar - DfromB",
1325  kVarForSparseFD,nbinsFD,xminFD,xmaxFD);
1326  fHistMassPtImpPar[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1327  "Mass vs. pt vs.true imppar -DfromB",
1328  kVarForSparseFD,nbinsFD,xminFD,xmaxFD);
1329  fHistMassPtImpPar[4]=new THnSparseF("hMassPtImpParBkg",
1330  "Mass vs. pt vs.imppar - backgr.",
1331  kVarForSparse,nbins,xmin,xmax);
1332 
1333 
1334  for(Int_t i=0; i<5; i++){
1335  for(Int_t iax=0; iax<kVarForSparse; iax++) fHistMassPtImpPar[i]->GetAxis(iax)->SetTitle(axTit[iax].Data());
1336  if(i == 2 || i == 3) fHistMassPtImpPar[i]->GetAxis(kVarForSparseFD-1)->SetTitle("p_{T}^{B} (GeV/c)");
1337  fOutput->Add(fHistMassPtImpPar[i]);
1338  }
1339 }
1340 
1341 //________________________________________________________________________
1344 
1345  const Int_t nVarPrompt = 2;
1346  const Int_t nVarFD = 3;
1347 
1348  Int_t nbinsPrompt[nVarPrompt]={200,100};
1349  Int_t nbinsFD[nVarFD]={200,100,200};
1350 
1351  Double_t xminPrompt[nVarPrompt] = {0.,-1.};
1352  Double_t xmaxPrompt[nVarPrompt] = {40.,1.};
1353 
1354  Double_t xminFD[nVarFD] = {0.,-1.,0.};
1355  Double_t xmaxFD[nVarFD] = {40.,1.,40.};
1356 
1357  //pt, y
1358  fMCAccPrompt = new THnSparseF("hMCAccPrompt","kStepMCAcceptance pt vs. y - promptD",nVarPrompt,nbinsPrompt,xminPrompt,xmaxPrompt);
1359  fMCAccPrompt->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
1360  fMCAccPrompt->GetAxis(1)->SetTitle("y");
1361 
1362  //pt,y,ptB
1363  fMCAccBFeed = new THnSparseF("hMCAccBFeed","kStepMCAcceptance pt vs. y vs. ptB - DfromB",nVarFD,nbinsFD,xminFD,xmaxFD);
1364  fMCAccBFeed->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
1365  fMCAccBFeed->GetAxis(1)->SetTitle("y");
1366  fMCAccBFeed->GetAxis(2)->SetTitle("p_{T}^{B} (GeV/c)");
1367 
1368  fOutput->Add(fMCAccPrompt);
1369  fOutput->Add(fMCAccBFeed);
1370 }
1371 
1372 //________________________________________________________________________
1373 void AliAnalysisTaskSEDplus::FillMCAcceptanceHistos(TClonesArray *arrayMC, AliAODMCHeader *mcHeader){
1375 
1376  const Int_t nProng = 3;
1377 
1378  Double_t zMCVertex = mcHeader->GetVtxZ(); //vertex MC
1379 
1380  for(Int_t iPart=0; iPart<arrayMC->GetEntriesFast(); iPart++){
1381  AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(arrayMC->At(iPart));
1382  if (TMath::Abs(mcPart->GetPdgCode()) == 411){
1383 
1384  Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,mcPart,fUseQuarkTagInKine);//Prompt = 4, FeedDown = 5
1385 
1386  Int_t deca = 0;
1387  Bool_t isGoodDecay=kFALSE;
1388  Int_t labDau[4]={-1,-1,-1,-1};
1389  Bool_t isInAcc = kFALSE;
1390  Bool_t isFidAcc = kFALSE;
1391 
1392  deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,mcPart,labDau);
1393  if(deca > 0) isGoodDecay=kTRUE;
1394 
1395  if(labDau[0]==-1){
1396  continue; //protection against unfilled array of labels
1397  }
1398 
1399  isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(mcPart->Pt(),mcPart->Y());
1400  isInAcc=CheckAcc(arrayMC,nProng,labDau);
1401 
1402  if(isGoodDecay && TMath::Abs(zMCVertex) < fRDCutsAnalysis->GetMaxVtxZ() && isFidAcc && isInAcc) {
1403  //for prompt
1404  if(orig == 4){
1405  //fill histo for prompt
1406  Double_t arrayMCprompt[2] = {mcPart->Pt(),mcPart->Y()};
1407  fMCAccPrompt->Fill(arrayMCprompt);
1408  }
1409  //for FD
1410  else if(orig == 5){
1411  Double_t ptB = AliVertexingHFUtils::GetBeautyMotherPt(arrayMC,mcPart);
1412  //fill histo for FD
1413  Double_t arrayMCFD[3] = {mcPart->Pt(),mcPart->Y(),ptB};
1414  fMCAccBFeed->Fill(arrayMCFD);
1415  }
1416  else
1417  continue;
1418  }
1419  }
1420  }
1421 }
1422 
1423 //________________________________________________________________________
1424 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
1425 {
1427  //
1428  if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
1429 
1430  fOutput = dynamic_cast<TList*> (GetOutputData(1));
1431  if (!fOutput) {
1432  printf("ERROR: fOutput not available\n");
1433  return;
1434  }
1435 
1436  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
1437  if(fHistNEvents){
1438  printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
1439  }else{
1440  printf("ERROR: fHistNEvents not available\n");
1441  return;
1442  }
1443 
1444  return;
1445 }
1446 //_________________________________________________________________________________________________
1447 Float_t AliAnalysisTaskSEDplus::GetTrueImpactParameter(const AliAODMCHeader *mcHeader, TClonesArray* arrayMC, const AliAODMCParticle *partDp) const {
1449 
1450  Double_t vtxTrue[3];
1451  mcHeader->GetVertex(vtxTrue);
1452  Double_t origD[3];
1453  partDp->XvYvZv(origD);
1454  Short_t charge=partDp->Charge();
1455  Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
1456  for(Int_t iDau=0; iDau<3; iDau++){
1457  pXdauTrue[iDau]=0.;
1458  pYdauTrue[iDau]=0.;
1459  pZdauTrue[iDau]=0.;
1460  }
1461 
1462  Int_t nDau=partDp->GetNDaughters();
1463  Int_t labelFirstDau = partDp->GetDaughter(0);
1464  if(nDau==3){
1465  for(Int_t iDau=0; iDau<3; 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  pXdauTrue[iDau]=part->Px();
1473  pYdauTrue[iDau]=part->Py();
1474  pZdauTrue[iDau]=part->Pz();
1475  }
1476  }else if(nDau==2){
1477  Int_t theDau=0;
1478  for(Int_t iDau=0; iDau<2; iDau++){
1479  Int_t ind = labelFirstDau+iDau;
1480  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1481  if(!part){
1482  AliError("Daughter particle not found in MC array");
1483  return 99999.;
1484  }
1485  Int_t pdgCode=TMath::Abs(part->GetPdgCode());
1486  if(pdgCode==211 || pdgCode==321){
1487  pXdauTrue[theDau]=part->Px();
1488  pYdauTrue[theDau]=part->Py();
1489  pZdauTrue[theDau]=part->Pz();
1490  ++theDau;
1491  }else{
1492  Int_t nDauRes=part->GetNDaughters();
1493  if(nDauRes==2){
1494  Int_t labelFirstDauRes = part->GetDaughter(0);
1495  for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
1496  Int_t indDR = labelFirstDauRes+iDauRes;
1497  AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
1498  if(!partDR){
1499  AliError("Daughter particle not found in MC array");
1500  return 99999.;
1501  }
1502 
1503  Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
1504  if(pdgCodeDR==211 || pdgCodeDR==321){
1505  pXdauTrue[theDau]=partDR->Px();
1506  pYdauTrue[theDau]=partDR->Py();
1507  pZdauTrue[theDau]=partDR->Pz();
1508  ++theDau;
1509  }
1510  }
1511  }
1512  }
1513  }
1514  if(theDau!=3){
1515  AliError("Wrong number of decay prongs");
1516  return 99999.;
1517  }
1518  }
1519 
1520  Double_t d0dummy[3]={0.,0.,0.};
1521  AliAODRecoDecayHF aodDplusMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
1522  return aodDplusMC.ImpParXY();
1523 
1524 }
1525 //_________________________________________________________________________________________________
1526 Float_t AliAnalysisTaskSEDplus::GetStrangenessWeights(const AliAODRecoDecayHF3Prong* d, TClonesArray* arrayMC, Float_t factor[3]) const {
1528 
1529  for(Int_t iprong=0;iprong<3;iprong++){
1530  factor[iprong]=1;
1531  AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
1532  Int_t labd= trad->GetLabel();
1533  if(labd>=0){
1534  AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
1535  if(dau){
1536  Int_t labm = dau->GetMother();
1537  if(labm>=0){
1538  AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
1539  if(mot){
1540  if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
1541  if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
1542  else factor[iprong]=1./.6;
1543  // fNentries->Fill(12);
1544  }
1545  if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
1546  factor[iprong]=1./0.25;
1547  // fNentries->Fill(13);
1548  }//if 3122
1549  }//if(mot)
1550  }//if labm>0
1551  }//if(dau)
1552  }//if labd>=0
1553  }//prong loop
1554 
1555  Float_t fact=1.;
1556  for(Int_t k=0;k<3;k++)fact=fact*factor[k];
1557  return fact;
1558 
1559 }
1560 
1561 //_________________________________________________________________
1562 Bool_t AliAnalysisTaskSEDplus::CheckAcc(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
1564  for (Int_t iProng = 0; iProng<nProng; iProng++){
1565  AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
1566  if(!mcPartDaughter) return kFALSE;
1567  Double_t eta = mcPartDaughter->Eta();
1568  Double_t pt = mcPartDaughter->Pt();
1569  if (TMath::Abs(eta) > 0.9 || pt < 0.1) return kFALSE;
1570  }
1571  return kTRUE;
1572 }
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
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:261
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: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)