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