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