AliPhysics  vAN-20150924 (e816f45)
 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(kFALSE),
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,Bool_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){
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  return;
749 }
750 
751 //________________________________________________________________________
752 void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
753 {
756 
757  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
758 
759  TClonesArray *array3Prong = 0;
760  TClonesArray *arrayLikeSign =0;
761  if(!aod && AODEvent() && IsStandardAOD()) {
762  // In case there is an AOD handler writing a standard AOD, use the AOD
763  // event in memory rather than the input (ESD) event.
764  aod = dynamic_cast<AliAODEvent*> (AODEvent());
765  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
766  // have to taken from the AOD event hold by the AliAODExtension
767  AliAODHandler* aodHandler = (AliAODHandler*)
768  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
769  if(aodHandler->GetExtensions()) {
770  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
771  AliAODEvent *aodFromExt = ext->GetAOD();
772  array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
773  arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
774  }
775  } else if(aod) {
776  array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
777  arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
778  }
779 
780  if(!aod || !array3Prong) {
781  printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
782  return;
783  }
784  if(!arrayLikeSign && fDoLS) {
785  printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
786  return;
787  }
788 
789  // fix for temporary bug in ESDfilter
790  // the AODs with null vertex pointer didn't pass the PhysSel
791  if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
792 
793  //Store the event in AliNormalizationCounter->To disable for Pb-Pb? Add a flag to disable it?
795 
796  fHistNEvents->Fill(0); // count event
797  Int_t runNumber=aod->GetRunNumber();
798 
799  //Event selection
800  Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
801  Float_t ntracks=aod->GetNumberOfTracks();
802  Float_t evCentr=0;
804  fHistCentrality[0]->Fill(ntracks,evCentr);
805  if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(2);
806  if(fRDCutsAnalysis->GetWhyRejection()==1) fHistNEvents->Fill(3);
807  if(fRDCutsAnalysis->GetWhyRejection()==2){fHistNEvents->Fill(4);fHistCentrality[2]->Fill(ntracks,evCentr);}
808  if(fRDCutsAnalysis->GetWhyRejection()==6)fHistNEvents->Fill(5);
809  if(fRDCutsAnalysis->GetWhyRejection()==7)fHistNEvents->Fill(6);
810 
811  TClonesArray *arrayMC=0;
812  AliAODMCHeader *mcHeader=0;
813  // load MC particles
814  if(fReadMC){
815  arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
816  if(!arrayMC) {
817  printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
818  return;
819  }
820 
821  // load MC header
822  mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
823  if(!mcHeader) {
824  printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
825  return;
826  }
827  }
828  if(fReadMC && fStepMCAcc){
829  if(aod->GetTriggerMask()==0 &&
830  (runNumber>=195344 && runNumber<=195677)){
831  // protection for events with empty trigger mask in p-Pb
832  return;
833  }
835  // events not passing the centrality selection can be removed immediately.
836 
837  FillMCAcceptanceHistos(arrayMC, mcHeader);
838  }
839  // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
840  //TString trigclass=aod->GetFiredTriggerClasses();
841  // Post the data already here
842  PostData(1,fOutput);
843  if(!isEvSel)return;
844  Int_t tracklets=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
845  // printf("ntracklet===%d\n",tracklets);
846  fSPDMult->Fill(tracklets);
847 
848  fHistCentrality[1]->Fill(ntracks,evCentr);
849  fHistNEvents->Fill(1);
850 
851  // AOD primary vertex
852  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
853  // vtx1->Print();
854  // TString primTitle = vtx1->GetTitle();
855  //if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0)fHistNEvents->Fill(2);
856 
857 
858  Int_t n3Prong = array3Prong->GetEntriesFast();
859  // printf("Number of D+->Kpipi: %d and of tracks: %d\n",n3Prong,aod->GetNumberOfTracks());
860 
861  Int_t nOS=0;
862  Int_t index;
863  Int_t pdgDgDplustoKpipi[3]={321,211,211};
864 
865  if(fDoLS>1){//Normalizations for LS
866  for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
867  AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
870  }
871  }
872  }else{//Standard analysis
873  // 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
874  //Double_t *cutsDplus = new (Double_t*)fRDCuts->GetCuts();
875  Int_t nSelectedloose=0,nSelectedtight=0;
876  for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
877  AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
878  fHistNEvents->Fill(7);
880  fHistNEvents->Fill(8);
881  continue;
882  }
883 
884  Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
885 
886  if(!fRDCutsAnalysis->GetIsSelectedCuts()) continue;
887 
888  Double_t etaD=d->Eta();
889  Double_t phiD=d->Phi();
890  if(fEtaSelection!=0){
891  if(fEtaSelection==1 && etaD<0) continue;
892  if(fEtaSelection==-1 && etaD>0) continue;
893  }
894 
895  Bool_t unsetvtx=kFALSE;
896  if(!d->GetOwnPrimaryVtx()){
897  d->SetOwnPrimaryVtx(vtx1);
898  unsetvtx=kTRUE;
899  }
900 
901  Double_t ptCand = d->Pt();
902  Int_t iPtBin = fRDCutsAnalysis->PtBin(ptCand);
903 
904  Bool_t recVtx=kFALSE;
905  AliAODVertex *origownvtx=0x0;
907  if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
908  if(fRDCutsAnalysis->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
909  else fRDCutsAnalysis->CleanOwnPrimaryVtx(d,aod,origownvtx);
910  }
911 
912  Int_t labDp=-1;
913  Bool_t isPrimary=kFALSE;
914  Bool_t isFeeddown=kFALSE;
915  Float_t pdgCode=-2;
916  Float_t trueImpParXY=0.;
917  Double_t ptB=-1.5;
918  if(fReadMC){
919  labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
920  if(labDp>=0){
921  AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
922  Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,partDp,fUseQuarkTagInKine);//Prompt = 4, FeedDown = 5
923  pdgCode=TMath::Abs(partDp->GetPdgCode());
924  if(orig==4){
925  isPrimary=kTRUE;
926  isFeeddown=kFALSE;
927  }else if(orig==5){
928  isPrimary=kFALSE;
929  isFeeddown=kTRUE;
930  trueImpParXY=GetTrueImpactParameter(mcHeader,arrayMC,partDp)*10000.;
931  ptB=AliVertexingHFUtils::GetBeautyMotherPt(arrayMC,partDp);
932  }else{
933  pdgCode=-3;
934  }
935  }else{
936  pdgCode=-1;
937  }
938  }
939 
940  Double_t invMass=d->InvMassDplus();
941  Double_t rapid=d->YDplus();
942  fYVsPt->Fill(ptCand,rapid,invMass);
943  if(passTightCuts) {fYVsPtTC->Fill(ptCand,rapid,invMass);nOS++;}
944  Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
945  if(isFidAcc){
946  fPtVsMass->Fill(invMass,ptCand);
947  if(passTightCuts){
948  fPtVsMassTC->Fill(invMass,ptCand);
949  fPhiEtaCand->Fill(etaD,phiD);
950  if(TMath::Abs(invMass-1.8696)<0.05) fPhiEtaCandSigReg->Fill(etaD,phiD);
951  }
952  }
953 
954 
955  Double_t dlen=0,cosp=0,maxdca=0,sigvert=0,sumD02=0,ptmax=0,dlxy=0,cxy=0, dxy=0;
957  dlen=d->DecayLength();
958  cosp=d->CosPointingAngle();
959  sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
960  maxdca=-9999.;
961  for(Int_t idau=0;idau<3;idau++) if(d->GetDCA(idau)>maxdca) maxdca=d->GetDCA(idau);
962  sigvert=d->GetSigmaVert();
963  ptmax=0;
964  for(Int_t i=0;i<3;i++){
965  if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
966  }
967  dlxy=d->NormalizedDecayLengthXY();
968  dxy = d->DecayLengthXY();
969  cxy=d->CosPointingAngleXY();
970  }
971  Double_t impparXY=d->ImpParXY()*10000.;
972  //for all THnSParse except for FD
973  Double_t arrayForSparse[6]={invMass,ptCand,impparXY,cxy,dxy,dlxy};
974 
975  //for THnSparse for FD
976  Double_t arrayForSparseFD[7]={invMass,ptCand,impparXY,cxy,dxy,dlxy,ptB};
977  Double_t arrayForSparseTrue[7]={invMass,ptCand,trueImpParXY,cxy,dxy,dlxy,ptB};
978 
979  //Ntuple
980  Float_t tmp[31];
981  if(fFillNtuple){
982  tmp[0]=pdgCode;
983  if(isFeeddown) tmp[0]+=5000.;
984  tmp[1]=d->Px();
985  tmp[2]=d->Py();
986  tmp[3]=d->Pz();
987  tmp[4]=d->Pt();
988  tmp[5]=d->GetCharge();
989  // tmp[5]=fRDCutsAnalysis->GetPIDBitMask(d);
990  tmp[6]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(0));
991  tmp[7]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(1));
992  tmp[8]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(2));
993  tmp[9]=d->PtProng(0);
994  tmp[10]=d->PtProng(1);
995  tmp[11]=d->PtProng(2);
996  tmp[12]=d->PProng(0);
997  tmp[13]=d->PProng(1);
998  tmp[14]=d->PProng(2);
999  tmp[15]=cosp;
1000  tmp[16]=cxy;
1001  tmp[17]=dlen;
1002  tmp[18]=d->NormalizedDecayLength();
1003  tmp[19]=d->DecayLengthXY();
1004  tmp[20]=dlxy;
1005  tmp[21]=d->InvMassDplus();
1006  tmp[22]=sigvert;
1007  tmp[23]=d->Getd0Prong(0);
1008  tmp[24]=d->Getd0Prong(1);
1009  tmp[25]=d->Getd0Prong(2);
1010  tmp[26]=maxdca;
1011  tmp[27]=ntracks;
1012  tmp[28]=fRDCutsAnalysis->GetCentrality(aod);
1013  tmp[29]=runNumber;
1014  tmp[30]=d->HasBadDaughters();
1015  fNtupleDplus->Fill(tmp);
1016  PostData(4,fNtupleDplus);
1017  }
1018 
1019  //Fill histos
1020  index=GetHistoIndex(iPtBin);
1021  if(isFidAcc){
1022  fHistNEvents->Fill(9);
1023  nSelectedloose++;
1024  fMassHist[index]->Fill(invMass);
1025  if(fCutsDistr){
1026  fCosPHist[index]->Fill(cosp);
1027  fDLenHist[index]->Fill(dlen);
1028  fSumd02Hist[index]->Fill(sumD02);
1029  fSigVertHist[index]->Fill(sigvert);
1030  fPtMaxHist[index]->Fill(ptmax);
1031  fPtKHist[index]->Fill(d->PtProng(1));
1032  fPtpi1Hist[index]->Fill(d->PtProng(0));
1033  fPtpi2Hist[index]->Fill(d->PtProng(2));
1034  fDCAHist[index]->Fill(maxdca);
1035  fDLxy[index]->Fill(dlxy);
1036  fCosxy[index]->Fill(cxy);
1037  fCorreld0Kd0pi[0]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
1038  d->Getd0Prong(2)*d->Getd0Prong(1));
1039  }
1040  if(passTightCuts){
1041  fHistNEvents->Fill(10);
1042  nSelectedtight++;
1043  fMassHistTC[index]->Fill(invMass);
1044  if(fCutsDistr){
1045  fDLxyTC[index]->Fill(dlxy);
1046  fCosxyTC[index]->Fill(cxy);
1047  }
1048  if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
1049  else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
1050  if(fDoImpPar){
1051  fHistMassPtImpParTC[0]->Fill(arrayForSparse);
1052  }
1053  }
1054  }
1055 
1056  if(fReadMC){
1057  if(isFidAcc){
1058  if(labDp>=0) {
1059  index=GetSignalHistoIndex(iPtBin);
1060  if(passTightCuts&&fDoImpPar){
1061  if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse);
1062  else if(isFeeddown){
1063  fHistMassPtImpParTC[2]->Fill(arrayForSparseFD);
1064  fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue);
1065  }
1066  }
1067  }else{
1068  index=GetBackgroundHistoIndex(iPtBin);
1069  if(passTightCuts&&fDoImpPar)fHistMassPtImpParTC[4]->Fill(arrayForSparse);
1070  }
1071 
1072  fMassHist[index]->Fill(invMass);
1073  if(fCutsDistr){
1074  Float_t fact=1.;
1075  Float_t factor[3]={1.,1.,1.};
1076  if(fUseStrangeness) fact=GetStrangenessWeights(d,arrayMC,factor);
1077  fCosPHist[index]->Fill(cosp,fact);
1078  fDLenHist[index]->Fill(dlen,fact);
1079  fDLxy[index]->Fill(dlxy);
1080  fCosxy[index]->Fill(cxy);
1081 
1082  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];
1083  fSumd02Hist[index]->Fill(sumd02s);
1084  fSigVertHist[index]->Fill(sigvert,fact);
1085  fPtMaxHist[index]->Fill(ptmax,fact);
1086  fPtKHist[index]->Fill(d->PtProng(1),fact);
1087  fPtpi1Hist[index]->Fill(d->PtProng(0),fact);
1088  fPtpi2Hist[index]->Fill(d->PtProng(2),fact);
1089  fDCAHist[index]->Fill(maxdca,fact);
1090  fCorreld0Kd0pi[1]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
1091  d->Getd0Prong(2)*d->Getd0Prong(1));
1092  }
1093  if(passTightCuts){
1094  fMassHistTC[index]->Fill(invMass);
1095  if(fCutsDistr){
1096  fDLxyTC[index]->Fill(dlxy);
1097  fCosxyTC[index]->Fill(cxy);
1098  }
1099  if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
1100  else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
1101  }
1102  }else{//outside fidAcc
1103  if(labDp>=0){
1104  fYVsPtSig->Fill(ptCand,rapid, invMass);
1105  if(passTightCuts)fYVsPtSigTC->Fill(ptCand,rapid, invMass);
1106  }
1107  }
1108  }//readmc
1109 
1110  if(recVtx)fRDCutsAnalysis->CleanOwnPrimaryVtx(d,aod,origownvtx);
1111 
1112  if(unsetvtx) d->UnsetOwnPrimaryVtx();
1113  }
1114  fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
1115  fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
1116  }
1117  //start LS analysis
1118  if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
1119 
1120  PostData(1,fOutput);
1121  PostData(2,fListCuts);
1122  PostData(3,fCounter);
1123  return;
1124 }
1125 
1126 //________________________________________________________________________
1129 
1130  TString hisname;
1131  Int_t indexLS=0;
1132  Int_t index=0;
1133  Int_t nbins=GetNBinsHistos();
1134  for(Int_t i=0;i<fNPtBins;i++){
1135 
1136  index=GetHistoIndex(i);
1137  indexLS=GetLSHistoIndex(i);
1138 
1139  hisname.Form("hLSPt%dLC",i);
1140  fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1141  fMassHistLS[indexLS]->Sumw2();
1142  hisname.Form("hLSPt%dTC",i);
1143  fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1144  fMassHistLSTC[indexLS]->Sumw2();
1145 
1146  hisname.Form("hCosPAllPt%dLS",i);
1147  fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1148  fCosPHistLS[index]->Sumw2();
1149  hisname.Form("hDLenAllPt%dLS",i);
1150  fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1151  fDLenHistLS[index]->Sumw2();
1152  hisname.Form("hSumd02AllPt%dLS",i);
1153  fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1154  fSumd02HistLS[index]->Sumw2();
1155  hisname.Form("hSigVertAllPt%dLS",i);
1156  fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1157  fSigVertHistLS[index]->Sumw2();
1158  hisname.Form("hPtMaxAllPt%dLS",i);
1159  fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1160  fPtMaxHistLS[index]->Sumw2();
1161  hisname.Form("hDCAAllPt%dLS",i);
1162  fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1163  fDCAHistLS[index]->Sumw2();
1164 
1165  index=GetSignalHistoIndex(i);
1166  indexLS++;
1167 
1168  hisname.Form("hLSPt%dLCnw",i);
1169  fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1170  fMassHistLS[indexLS]->Sumw2();
1171  hisname.Form("hLSPt%dTCnw",i);
1172  fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1173  fMassHistLSTC[indexLS]->Sumw2();
1174 
1175  hisname.Form("hCosPSigPt%dLS",i);
1176  fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1177  fCosPHistLS[index]->Sumw2();
1178  hisname.Form("hDLenSigPt%dLS",i);
1179  fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1180  fDLenHistLS[index]->Sumw2();
1181  hisname.Form("hSumd02SigPt%dLS",i);
1182  fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1183  fSumd02HistLS[index]->Sumw2();
1184  hisname.Form("hSigVertSigPt%dLS",i);
1185  fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1186  fSigVertHistLS[index]->Sumw2();
1187  hisname.Form("hPtMaxSigPt%dLS",i);
1188  fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1189  fPtMaxHistLS[index]->Sumw2();
1190  hisname.Form("hDCASigPt%dLS",i);
1191  fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1192  fDCAHistLS[index]->Sumw2();
1193 
1194  index=GetBackgroundHistoIndex(i);
1195  indexLS++;
1196 
1197  hisname.Form("hLSPt%dLCntrip",i);
1198  fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1199  fMassHistLS[indexLS]->Sumw2();
1200  hisname.Form("hLSPt%dTCntrip",i);
1201  fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1202  fMassHistLSTC[indexLS]->Sumw2();
1203 
1204  hisname.Form("hCosPBkgPt%dLS",i);
1205  fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1206  fCosPHistLS[index]->Sumw2();
1207  hisname.Form("hDLenBkgPt%dLS",i);
1208  fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1209  fDLenHistLS[index]->Sumw2();
1210  hisname.Form("hSumd02BkgPt%dLS",i);
1211  fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1212  fSumd02HistLS[index]->Sumw2();
1213  hisname.Form("hSigVertBkgPt%dLS",i);
1214  fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1215  fSigVertHistLS[index]->Sumw2();
1216  hisname.Form("hPtMaxBkgPt%dLS",i);
1217  fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1218  fPtMaxHistLS[index]->Sumw2();
1219  hisname.Form("hDCABkgPt%dLS",i);
1220  fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1221  fDCAHistLS[index]->Sumw2();
1222 
1223  indexLS++;
1224  hisname.Form("hLSPt%dLCntripsinglecut",i);
1225  fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1226  fMassHistLS[indexLS]->Sumw2();
1227  hisname.Form("hLSPt%dTCntripsinglecut",i);
1228  fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1229  fMassHistLSTC[indexLS]->Sumw2();
1230 
1231  indexLS++;
1232  hisname.Form("hLSPt%dLCspc",i);
1233  fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1234  fMassHistLS[indexLS]->Sumw2();
1235  hisname.Form("hLSPt%dTCspc",i);
1236  fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1237  fMassHistLSTC[indexLS]->Sumw2();
1238  }
1239 
1240  for(Int_t i=0; i<3*fNPtBins; i++){
1241  fOutput->Add(fCosPHistLS[i]);
1242  fOutput->Add(fDLenHistLS[i]);
1243  fOutput->Add(fSumd02HistLS[i]);
1244  fOutput->Add(fSigVertHistLS[i]);
1245  fOutput->Add(fPtMaxHistLS[i]);
1246  fOutput->Add(fDCAHistLS[i]);
1247  }
1248  for(Int_t i=0; i<5*fNPtBins; i++){
1249  fOutput->Add(fMassHistLS[i]);
1250  fOutput->Add(fMassHistLSTC[i]);
1251  }
1252 }
1253 
1254 //________________________________________________________________________
1257 
1258  Int_t nmassbins=GetNBinsHistos();
1259  Double_t maxmult;
1260  if(fSystem==1) maxmult=5000;
1261  else maxmult=200;
1262 
1263  const Int_t Nvar = 6;
1264  const Int_t NvarFD = 7;
1265 
1266  //dimensions for THnSparse which are NOT for BFeed
1267  Int_t nbins[Nvar]={nmassbins,80,fNImpParBins,30,100,30};
1268  Double_t xmin[Nvar]={fLowmasslimit,0.,fLowerImpPar,0.97,0.,0.};
1269  Double_t xmax[Nvar]={fUpmasslimit,40.,fHigherImpPar,1.,1,30.};
1270 
1271  //dimensions for THnSparse for BFeed
1272  Int_t nbinsFD[NvarFD]={nmassbins,80,fNImpParBins,30,100,30,84};
1273  Double_t xminFD[NvarFD]={fLowmasslimit,0.,fLowerImpPar,0.97,0.,0.,-2};
1274  Double_t xmaxFD[NvarFD]={fUpmasslimit,40.,fHigherImpPar,1.,1,30.,40};
1275 
1276  //mass, pt, imppar, cosPoinXY, decLXY, norm decLXY (for BFeed also ptB)
1277  //mass, pt, imppar, cosPoinXY, decLXY, norm decLXY (for BFeed also ptB)
1278  fHistMassPtImpParTC[0]=new THnSparseF("hMassPtImpParAll",
1279  "Mass vs. pt vs.imppar - All",
1280  Nvar,nbins,xmin,xmax);
1281  fHistMassPtImpParTC[1]=new THnSparseF("hMassPtImpParPrompt",
1282  "Mass vs. pt vs.imppar - promptD",
1283  Nvar,nbins,xmin,xmax);
1284  fHistMassPtImpParTC[2]=new THnSparseF("hMassPtImpParBfeed",
1285  "Mass vs. pt vs.imppar - DfromB",
1286  NvarFD,nbinsFD,xminFD,xmaxFD);
1287  fHistMassPtImpParTC[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1288  "Mass vs. pt vs.true imppar -DfromB",
1289  NvarFD,nbinsFD,xminFD,xmaxFD);
1290  fHistMassPtImpParTC[4]=new THnSparseF("hMassPtImpParBkg",
1291  "Mass vs. pt vs.imppar - backgr.",
1292  Nvar,nbins,xmin,xmax);
1293 
1294  for(Int_t i=0; i<5; i++){
1295  fHistMassPtImpParTC[i]->GetAxis(0)->SetTitle("M_{K#pi#pi} (GeV/c^{2})");
1296  fHistMassPtImpParTC[i]->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
1297  fHistMassPtImpParTC[i]->GetAxis(2)->SetTitle("Imp Par (#mum)");
1298  fHistMassPtImpParTC[i]->GetAxis(3)->SetTitle("cos(#theta_{P}^{xy})");
1299  fHistMassPtImpParTC[i]->GetAxis(4)->SetTitle("decL XY (cm)");
1300  fHistMassPtImpParTC[i]->GetAxis(5)->SetTitle("Norm decL XY");
1301 
1302  if(i == 2 || i == 3)
1303  fHistMassPtImpParTC[i]->GetAxis(6)->SetTitle("p_{T}^{B} (GeV/c)");
1304 
1305  fOutput->Add(fHistMassPtImpParTC[i]);
1306  }
1307 }
1308 
1309 //________________________________________________________________________
1312 
1313  const Int_t NvarPrompt = 2;
1314  const Int_t NvarFD = 3;
1315 
1316  Int_t nbinsPrompt[NvarPrompt]={200,100};
1317  Int_t nbinsFD[NvarFD]={200,100,200};
1318 
1319  Double_t xminPrompt[NvarPrompt] = {0.,-1.};
1320  Double_t xmaxPrompt[NvarPrompt] = {40.,1.};
1321 
1322  Double_t xminFD[NvarFD] = {0.,-1.,0.};
1323  Double_t xmaxFD[NvarFD] = {40.,1.,40.};
1324 
1325  //pt, y
1326  fMCAccPrompt = new THnSparseF("hMCAccPrompt","kStepMCAcceptance pt vs. y - promptD",NvarPrompt,nbinsPrompt,xminPrompt,xmaxPrompt);
1327  fMCAccPrompt->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
1328  fMCAccPrompt->GetAxis(1)->SetTitle("y");
1329 
1330  //pt,y,ptB
1331  fMCAccBFeed = new THnSparseF("hMCAccBFeed","kStepMCAcceptance pt vs. y vs. ptB - DfromB",NvarFD,nbinsFD,xminFD,xmaxFD);
1332  fMCAccBFeed->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
1333  fMCAccBFeed->GetAxis(1)->SetTitle("y");
1334  fMCAccBFeed->GetAxis(2)->SetTitle("p_{T}^{B} (GeV/c)");
1335 
1336  fOutput->Add(fMCAccPrompt);
1337  fOutput->Add(fMCAccBFeed);
1338 }
1339 
1340 //________________________________________________________________________
1341 void AliAnalysisTaskSEDplus::FillMCAcceptanceHistos(TClonesArray *arrayMC, AliAODMCHeader *mcHeader){
1343 
1344  const Int_t nProng = 3;
1345 
1346  Double_t zMCVertex = mcHeader->GetVtxZ(); //vertex MC
1347 
1348  for(Int_t iPart=0; iPart<arrayMC->GetEntriesFast(); iPart++){
1349  AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(arrayMC->At(iPart));
1350  if (TMath::Abs(mcPart->GetPdgCode()) == 411){
1351 
1352  Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,mcPart,fUseQuarkTagInKine);//Prompt = 4, FeedDown = 5
1353 
1354  Int_t deca = 0;
1355  Bool_t isGoodDecay=kFALSE;
1356  Int_t labDau[4]={-1,-1,-1,-1};
1357  Bool_t isInAcc = kFALSE;
1358  Bool_t isFidAcc = kFALSE;
1359 
1360  deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,mcPart,labDau);
1361  if(deca > 0) isGoodDecay=kTRUE;
1362 
1363  if(labDau[0]==-1){
1364  continue; //protection against unfilled array of labels
1365  }
1366 
1367  isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(mcPart->Pt(),mcPart->Y());
1368  isInAcc=CheckAcc(arrayMC,nProng,labDau);
1369 
1370  if(isGoodDecay && TMath::Abs(zMCVertex) < fRDCutsAnalysis->GetMaxVtxZ() && isFidAcc && isInAcc) {
1371  //for prompt
1372  if(orig == 4){
1373  //fill histo for prompt
1374  Double_t arrayMCprompt[2] = {mcPart->Pt(),mcPart->Y()};
1375  fMCAccPrompt->Fill(arrayMCprompt);
1376  }
1377  //for FD
1378  else if(orig == 5){
1379  Double_t ptB = AliVertexingHFUtils::GetBeautyMotherPt(arrayMC,mcPart);
1380  //fill histo for FD
1381  Double_t arrayMCFD[3] = {mcPart->Pt(),mcPart->Y(),ptB};
1382  fMCAccBFeed->Fill(arrayMCFD);
1383  }
1384  else
1385  continue;
1386  }
1387  }
1388  }
1389 }
1390 
1391 //________________________________________________________________________
1392 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
1393 {
1395  //
1396  if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
1397 
1398  fOutput = dynamic_cast<TList*> (GetOutputData(1));
1399  if (!fOutput) {
1400  printf("ERROR: fOutput not available\n");
1401  return;
1402  }
1403 
1404  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
1405  if(fHistNEvents){
1406  printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
1407  }else{
1408  printf("ERROR: fHistNEvents not available\n");
1409  return;
1410  }
1411 
1412  return;
1413 }
1414 //_________________________________________________________________________________________________
1415 Float_t AliAnalysisTaskSEDplus::GetTrueImpactParameter(const AliAODMCHeader *mcHeader, TClonesArray* arrayMC, const AliAODMCParticle *partDp) const {
1417 
1418  Double_t vtxTrue[3];
1419  mcHeader->GetVertex(vtxTrue);
1420  Double_t origD[3];
1421  partDp->XvYvZv(origD);
1422  Short_t charge=partDp->Charge();
1423  Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
1424  for(Int_t iDau=0; iDau<3; iDau++){
1425  pXdauTrue[iDau]=0.;
1426  pYdauTrue[iDau]=0.;
1427  pZdauTrue[iDau]=0.;
1428  }
1429 
1430  Int_t nDau=partDp->GetNDaughters();
1431  Int_t labelFirstDau = partDp->GetDaughter(0);
1432  if(nDau==3){
1433  for(Int_t iDau=0; iDau<3; iDau++){
1434  Int_t ind = labelFirstDau+iDau;
1435  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1436  if(!part){
1437  AliError("Daughter particle not found in MC array");
1438  return 99999.;
1439  }
1440  pXdauTrue[iDau]=part->Px();
1441  pYdauTrue[iDau]=part->Py();
1442  pZdauTrue[iDau]=part->Pz();
1443  }
1444  }else if(nDau==2){
1445  Int_t theDau=0;
1446  for(Int_t iDau=0; iDau<2; iDau++){
1447  Int_t ind = labelFirstDau+iDau;
1448  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1449  if(!part){
1450  AliError("Daughter particle not found in MC array");
1451  return 99999.;
1452  }
1453  Int_t pdgCode=TMath::Abs(part->GetPdgCode());
1454  if(pdgCode==211 || pdgCode==321){
1455  pXdauTrue[theDau]=part->Px();
1456  pYdauTrue[theDau]=part->Py();
1457  pZdauTrue[theDau]=part->Pz();
1458  ++theDau;
1459  }else{
1460  Int_t nDauRes=part->GetNDaughters();
1461  if(nDauRes==2){
1462  Int_t labelFirstDauRes = part->GetDaughter(0);
1463  for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
1464  Int_t indDR = labelFirstDauRes+iDauRes;
1465  AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
1466  if(!partDR){
1467  AliError("Daughter particle not found in MC array");
1468  return 99999.;
1469  }
1470 
1471  Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
1472  if(pdgCodeDR==211 || pdgCodeDR==321){
1473  pXdauTrue[theDau]=partDR->Px();
1474  pYdauTrue[theDau]=partDR->Py();
1475  pZdauTrue[theDau]=partDR->Pz();
1476  ++theDau;
1477  }
1478  }
1479  }
1480  }
1481  }
1482  if(theDau!=3){
1483  AliError("Wrong number of decay prongs");
1484  return 99999.;
1485  }
1486  }
1487 
1488  Double_t d0dummy[3]={0.,0.,0.};
1489  AliAODRecoDecayHF aodDplusMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
1490  return aodDplusMC.ImpParXY();
1491 
1492 }
1493 //_________________________________________________________________________________________________
1494 Float_t AliAnalysisTaskSEDplus::GetStrangenessWeights(const AliAODRecoDecayHF3Prong* d, TClonesArray* arrayMC, Float_t factor[3]) const {
1496 
1497  for(Int_t iprong=0;iprong<3;iprong++){
1498  factor[iprong]=1;
1499  AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
1500  Int_t labd= trad->GetLabel();
1501  if(labd>=0){
1502  AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
1503  if(dau){
1504  Int_t labm = dau->GetMother();
1505  if(labm>=0){
1506  AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
1507  if(mot){
1508  if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
1509  if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
1510  else factor[iprong]=1./.6;
1511  // fNentries->Fill(12);
1512  }
1513  if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
1514  factor[iprong]=1./0.25;
1515  // fNentries->Fill(13);
1516  }//if 3122
1517  }//if(mot)
1518  }//if labm>0
1519  }//if(dau)
1520  }//if labd>=0
1521  }//prong loop
1522 
1523  Float_t fact=1.;
1524  for(Int_t k=0;k<3;k++)fact=fact*factor[k];
1525  return fact;
1526 
1527 }
1528 
1529 //_________________________________________________________________
1530 Bool_t AliAnalysisTaskSEDplus::CheckAcc(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
1532  for (Int_t iProng = 0; iProng<nProng; iProng++){
1533  AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
1534  if(!mcPartDaughter) return kFALSE;
1535  Double_t eta = mcPartDaughter->Eta();
1536  Double_t pt = mcPartDaughter->Pt();
1537  if (TMath::Abs(eta) > 0.9 || pt < 0.1) return kFALSE;
1538  }
1539  return kTRUE;
1540 }
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
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)
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)
Bool_t fFillNtuple
limits for the Pt bins
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)