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