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