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