]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added Init() for AliNormalizationCounter
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDplus.cxx
CommitLineData
d486095a 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
27de2dfb 16/* $Id$ */
17
993bfba1 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)
4afc48a2 23// Authors: Renu Bala, bala@to.infn.it
24// F. Prino, prino@to.infn.it
25// G. Ortona, ortona@to.infn.it
d486095a 26/////////////////////////////////////////////////////////////
27
28#include <TClonesArray.h>
29#include <TNtuple.h>
13808a30 30#include <TCanvas.h>
d486095a 31#include <TList.h>
1f4e9722 32#include <TString.h>
4afc48a2 33#include <TDatabasePDG.h>
b557eb43 34
35#include "AliAnalysisManager.h"
4c230f6d 36#include "AliRDHFCutsDplustoKpipi.h"
b557eb43 37#include "AliAODHandler.h"
d486095a 38#include "AliAODEvent.h"
39#include "AliAODVertex.h"
40#include "AliAODTrack.h"
d486095a 41#include "AliAODRecoDecayHF3Prong.h"
42#include "AliAnalysisVertexingHF.h"
43#include "AliAnalysisTaskSE.h"
44#include "AliAnalysisTaskSEDplus.h"
a96083b9 45#include "AliNormalizationCounter.h"
d486095a 46
47ClassImp(AliAnalysisTaskSEDplus)
48
49
50//________________________________________________________________________
51AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
52AliAnalysisTaskSE(),
595cc7e2 53 fOutput(0),
54 fHistNEvents(0),
55 fPtVsMass(0),
56 fPtVsMassTC(0),
57 fYVsPt(0),
58 fYVsPtTC(0),
59 fYVsPtSig(0),
60 fYVsPtSigTC(0),
61 fNtupleDplus(0),
62 fUpmasslimit(1.965),
63 fLowmasslimit(1.765),
64 fNPtBins(0),
65 fBinWidth(0.002),
66 fListCuts(0),
67 fRDCutsProduction(0),
68 fRDCutsAnalysis(0),
a96083b9 69 fCounter(0),
595cc7e2 70 fFillNtuple(kFALSE),
71 fReadMC(kFALSE),
a96083b9 72 fUseStrangeness(kFALSE),
5fc4893f 73 fUseBit(kTRUE),
2edd194d 74 fCutsDistr(kFALSE),
7d9ca2b5 75 fDoImpPar(kFALSE),
3de19caa 76 fNImpParBins(400),
77 fLowerImpPar(-2000.),
78 fHigherImpPar(2000.),
3ec9254b 79 fDoLS(0)
d486095a 80{
82bb8d4b 81 // Default constructor
d486095a 82}
83
84//________________________________________________________________________
4c230f6d 85AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,AliRDHFCutsDplustoKpipi *dpluscutsana,AliRDHFCutsDplustoKpipi *dpluscutsprod,Bool_t fillNtuple):
3de19caa 86 AliAnalysisTaskSE(name),
87 fOutput(0),
88 fHistNEvents(0),
89 fPtVsMass(0),
90 fPtVsMassTC(0),
91 fYVsPt(0),
92 fYVsPtTC(0),
93 fYVsPtSig(0),
94 fYVsPtSigTC(0),
95 fNtupleDplus(0),
96 fUpmasslimit(1.965),
97 fLowmasslimit(1.765),
98 fNPtBins(0),
99 fBinWidth(0.002),
100 fListCuts(0),
101 fRDCutsProduction(dpluscutsprod),
102 fRDCutsAnalysis(dpluscutsana),
103 fCounter(0),
104 fFillNtuple(fillNtuple),
105 fReadMC(kFALSE),
106 fUseStrangeness(kFALSE),
107 fUseBit(kTRUE),
2edd194d 108 fCutsDistr(kFALSE),
3de19caa 109 fDoImpPar(kFALSE),
110 fNImpParBins(400),
111 fLowerImpPar(-2000.),
112 fHigherImpPar(2000.),
113 fDoLS(0)
d486095a 114{
4c230f6d 115 //
116 // Standrd constructor
117 //
a6003e0a 118 fNPtBins=fRDCutsAnalysis->GetNPtBins();
d486095a 119 // Default constructor
82bb8d4b 120 // Output slot #1 writes into a TList container
d486095a 121 DefineOutput(1,TList::Class()); //My private output
4c230f6d 122 // Output slot #2 writes cut to private output
123 // DefineOutput(2,AliRDHFCutsDplustoKpipi::Class());
124 DefineOutput(2,TList::Class());
a96083b9 125// Output slot #3 writes cut to private output
126 DefineOutput(3,AliNormalizationCounter::Class());
127
1f4e9722 128 if(fFillNtuple){
a96083b9 129 // Output slot #4 writes into a TNtuple container
130 DefineOutput(4,TNtuple::Class()); //My private output
1f4e9722 131 }
d486095a 132}
133
134//________________________________________________________________________
135AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
136{
4c230f6d 137 //
d486095a 138 // Destructor
4c230f6d 139 //
d486095a 140 if (fOutput) {
141 delete fOutput;
142 fOutput = 0;
143 }
137c7bcd 144 if(fHistNEvents){
145 delete fHistNEvents;
146 fHistNEvents=0;
147 }
37fc80e6 148 for(Int_t i=0;i<3;i++){
149 if(fHistCentrality[i]){delete fHistCentrality[i]; fHistCentrality[i]=0;}
150 }
cdc197b1 151
137c7bcd 152 for(Int_t i=0;i<3*fNPtBins;i++){
153 if(fMassHist[i]){ delete fMassHist[i]; fMassHist[i]=0;}
154 if(fCosPHist[i]){ delete fCosPHist[i]; fCosPHist[i]=0;}
155 if(fDLenHist[i]){ delete fDLenHist[i]; fDLenHist[i]=0;}
156 if(fSumd02Hist[i]){ delete fSumd02Hist[i]; fSumd02Hist[i]=0;}
157 if(fSigVertHist[i]){ delete fSigVertHist[i]; fSigVertHist[i]=0;}
158 if(fPtMaxHist[i]){ delete fPtMaxHist[i]; fPtMaxHist[i]=0;}
3ec9254b 159 if(fPtKHist[i]){ delete fPtKHist[i]; fPtKHist[i]=0;}
160 if(fPtpi1Hist[i]){ delete fPtpi1Hist[i]; fPtpi1Hist[i]=0;}
161 if(fPtpi2Hist[i]){ delete fPtpi2Hist[i]; fPtpi2Hist[i]=0;}
137c7bcd 162 if(fDCAHist[i]){ delete fDCAHist[i]; fDCAHist[i]=0;}
163 if(fMassHistTC[i]){ delete fMassHistTC[i]; fMassHistTC[i]=0;}
bb69f127 164 if(fMassHistTCPlus[i]){ delete fMassHistTCPlus[i]; fMassHistTCPlus[i]=0;}
165 if(fMassHistTCMinus[i]){ delete fMassHistTCMinus[i]; fMassHistTCMinus[i]=0;}
4c230f6d 166
3ec9254b 167 if(fDLxy[i]){delete fDLxy[i]; fDLxy[i]=0;}
168 if(fDLxyTC[i]){delete fDLxyTC[i]; fDLxyTC[i]=0;}
169 if(fCosxy[i]){delete fCosxy[i]; fCosxy[i]=0;}
170 if(fCosxyTC[i]){delete fCosxyTC[i]; fCosxyTC[i]=0;}
137c7bcd 171 if(fMassHistLS[i]){ delete fMassHistLS[i]; fMassHistLS[i]=0;}
172 if(fCosPHistLS[i]){ delete fCosPHistLS[i]; fCosPHistLS[i]=0;}
173 if(fDLenHistLS[i]){ delete fDLenHistLS[i]; fDLenHistLS[i]=0;}
174 if(fSumd02HistLS[i]){ delete fSumd02HistLS[i]; fSumd02HistLS[i]=0;}
175 if(fSigVertHistLS[i]){ delete fSigVertHistLS[i]; fSigVertHistLS[i]=0;}
176 if(fPtMaxHistLS[i]){ delete fPtMaxHistLS[i]; fPtMaxHistLS[i]=0;}
177 if(fDCAHistLS[i]){ delete fDCAHistLS[i]; fDCAHistLS[i]=0;}
178 if(fMassHistLSTC[i]){ delete fMassHistLSTC[i]; fMassHistLSTC[i]=0;}
179 }
cdc197b1 180
181 for(Int_t i=0;i<3;i++){
182 if(fCorreld0Kd0pi[i]){ delete fCorreld0Kd0pi[i]; fCorreld0Kd0pi[i]=0;}
183 }
184
137c7bcd 185 if(fPtVsMass){
186 delete fPtVsMass;
187 fPtVsMass=0;
188 }
189 if(fPtVsMassTC){
190 delete fPtVsMassTC;
191 fPtVsMassTC=0;
192 }
193 if(fYVsPt){
194 delete fYVsPt;
195 fYVsPt=0;
196 }
197 if(fYVsPtTC){
198 delete fYVsPtTC;
199 fYVsPtTC=0;
200 }
201 if(fYVsPtSig){
202 delete fYVsPtSig;
203 fYVsPtSig=0;
204 }
205 if(fYVsPtSigTC){
206 delete fYVsPtSigTC;
207 fYVsPtSigTC=0;
208 }
209 if(fNtupleDplus){
210 delete fNtupleDplus;
211 fNtupleDplus=0;
212 }
4c230f6d 213 if (fListCuts) {
214 delete fListCuts;
215 fListCuts = 0;
d486095a 216 }
4c230f6d 217 if(fRDCutsProduction){
218 delete fRDCutsProduction;
219 fRDCutsProduction = 0;
220 }
137c7bcd 221 if(fRDCutsAnalysis){
4c230f6d 222 delete fRDCutsAnalysis;
223 fRDCutsAnalysis = 0;
224 }
7d9ca2b5 225 for(Int_t i=0; i<5; i++){
226 delete fHistMassPtImpParTC[i];
227 }
a96083b9 228 if(fCounter){
229 delete fCounter;
230 fCounter = 0;
231 }
137c7bcd 232
233
d486095a 234}
82bb8d4b 235//_________________________________________________________________
236void AliAnalysisTaskSEDplus::SetMassLimits(Float_t range){
4c230f6d 237 // set invariant mass limits
993bfba1 238 Float_t bw=GetBinWidth();
82bb8d4b 239 fUpmasslimit = 1.865+range;
240 fLowmasslimit = 1.865-range;
993bfba1 241 SetBinWidth(bw);
82bb8d4b 242}
243//_________________________________________________________________
244void AliAnalysisTaskSEDplus::SetMassLimits(Float_t lowlimit, Float_t uplimit){
4c230f6d 245 // set invariant mass limits
82bb8d4b 246 if(uplimit>lowlimit)
247 {
993bfba1 248 Float_t bw=GetBinWidth();
82bb8d4b 249 fUpmasslimit = lowlimit;
250 fLowmasslimit = uplimit;
993bfba1 251 SetBinWidth(bw);
82bb8d4b 252 }
253}
993bfba1 254//________________________________________________________________
255void AliAnalysisTaskSEDplus::SetBinWidth(Float_t w){
256 Float_t width=w;
257 Int_t nbins=(Int_t)((fUpmasslimit-fLowmasslimit)/width+0.5);
258 Int_t missingbins=4-nbins%4;
259 nbins=nbins+missingbins;
260 width=(fUpmasslimit-fLowmasslimit)/nbins;
261 if(missingbins!=0){
262 printf("AliAnalysisTaskSEDplus::SetBinWidth: W-bin width of %f will produce histograms not rebinnable by 4. New width set to %f\n",w,width);
263 }
264 else{
265 if(fDebug>1) printf("AliAnalysisTaskSEDplus::SetBinWidth: width set to %f\n",width);
266 }
267 fBinWidth=width;
268}
82bb8d4b 269//_________________________________________________________________
993bfba1 270Int_t AliAnalysisTaskSEDplus::GetNBinsHistos(){
271 return (Int_t)((fUpmasslimit-fLowmasslimit)/fBinWidth+0.5);
272}
82bb8d4b 273//_________________________________________________________________
274void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS){
4c230f6d 275 //
276 //
277 // Fill the Like Sign histograms
278 //
3ec9254b 279 if(fDebug>1)printf("started LS\n");
2edd194d 280
5fc4893f 281 //histograms for like sign
282 Int_t nbins=GetNBinsHistos();;
283 TH1F *histLSPlus = new TH1F("LSPlus","LSPlus",nbins,fLowmasslimit,fUpmasslimit);
284 TH1F *histLSMinus = new TH1F("LSMinus","LSMinus",nbins,fLowmasslimit,fUpmasslimit);
285
82bb8d4b 286 Int_t nPosTrks=0,nNegTrks=0;
82bb8d4b 287 Int_t nOStriplets = arrayOppositeSign->GetEntriesFast();
82bb8d4b 288 Int_t nDplusLS=0;
5fc4893f 289 Int_t nDminusLS=0;
82bb8d4b 290 Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
5fc4893f 291 Int_t index=0;
292
293 // loop over like sign candidates
82bb8d4b 294 for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
295 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
5fc4893f 296 if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts))continue;
82bb8d4b 297 Bool_t unsetvtx=kFALSE;
298 if(!d->GetOwnPrimaryVtx()) {
299 d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
300 unsetvtx=kTRUE;
301 }
4b99797c 302 Double_t ptCand = d->Pt();
303 Int_t iPtBin = fRDCutsAnalysis->PtBin(ptCand);
304 if(iPtBin<0)continue;
305 Int_t sign= d->GetCharge();
306 if(sign>0){
307 nPosTrks++;
308 }else{
309 nNegTrks++;
310 }
5fc4893f 311 // if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)){
312 fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
313 Int_t passTightCuts=fRDCutsAnalysis->GetIsSelectedCuts();
314 if(passTightCuts>0){
82bb8d4b 315 Float_t invMass = d->InvMassDplus();
82bb8d4b 316 index=GetLSHistoIndex(iPtBin);
82bb8d4b 317 fMassHistLS[index+1]->Fill(invMass);
5fc4893f 318 if(sign>0){
319 histLSPlus->Fill(invMass);
320 nDplusLS++;
321 }else{
322 histLSMinus->Fill(invMass);
323 nDminusLS++;
324 }
2edd194d 325 if(fCutsDistr){
5fc4893f 326 Double_t dlen=d->DecayLength();
327 Double_t cosp=d->CosPointingAngle();
328 Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
329 Double_t dca=d->GetDCA();
330 Double_t sigvert=d->GetSigmaVert();
331 Double_t ptmax=0;
332 for(Int_t i=0;i<3;i++){
333 if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
334 }
335 fCosPHistLS[iPtBin]->Fill(cosp);
336 fDLenHistLS[iPtBin]->Fill(dlen);
337 fSumd02HistLS[iPtBin]->Fill(sumD02);
338 fSigVertHistLS[iPtBin]->Fill(sigvert);
339 fPtMaxHistLS[iPtBin]->Fill(ptmax);
340 fDCAHistLS[iPtBin]->Fill(dca);
82bb8d4b 341 }
342 }
343 if(unsetvtx) d->UnsetOwnPrimaryVtx();
344 }
5fc4893f 345 //wei:normalized to the number of combinations (production)
346 //wei2:normalized to the number of LS/OS (production)
347 //wei3:normalized to the number of LS/OS (analysis)
348 //wei4:normalized to the number of combinations (analysis)
349 Float_t wei2=0;
350 if(nLikeSign!=0)wei2 = (Float_t)nOStriplets/(Float_t)nLikeSign;
351 Float_t wei3=0;
3ec9254b 352 if(nDplusLS!=0)wei3 = (Float_t)nDplusOS/(Float_t)(nDplusLS+nDminusLS);
5fc4893f 353 Float_t weiplus=1.,weiminus=1.;
354 Float_t wei4plus=1.,wei4minus=1.;
355 //wei* should be automatically protected, since to get a triplet there must be at least 3 good tracks in the event
356 if(nPosTrks>2)weiplus=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.);
357 if(nDplusLS>2)wei4plus=3.*(Float_t)nDminusLS/((Float_t)nDplusLS-2.);
358 if(nNegTrks>2)weiminus=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.);
359 if(nDminusLS>2)wei4minus=3.*(Float_t)nDplusLS/((Float_t)nDminusLS-2.);
360
361 fMassHistLS[index]->Add(histLSPlus,weiplus);
362 fMassHistLS[index]->Add(histLSMinus,weiminus);
363 fMassHistLS[index+2]->Add(histLSPlus,wei2);
364 fMassHistLS[index+2]->Add(histLSMinus,wei2);
365 fMassHistLS[index+3]->Add(histLSPlus,wei3);
366 fMassHistLS[index+3]->Add(histLSMinus,wei3);
367 fMassHistLS[index+4]->Add(histLSPlus,wei4plus);
368 fMassHistLS[index+4]->Add(histLSMinus,wei4minus);
369
370 delete histLSPlus;histLSPlus=0;
371 delete histLSMinus;histLSMinus=0;
82bb8d4b 372
5fc4893f 373 if(fDebug>1) printf("LS analysis terminated\n");
82bb8d4b 374}
d486095a 375
4c230f6d 376//__________________________________________
377void AliAnalysisTaskSEDplus::Init(){
378 //
379 // Initialization
380 //
d486095a 381 if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
4c230f6d 382
383 //PostData(2,fRDCutsloose);//we should then put those cuts in a tlist if we have more than 1
384 fListCuts=new TList();
90993728 385 AliRDHFCutsDplustoKpipi *production = new AliRDHFCutsDplustoKpipi(*fRDCutsProduction);
386 production->SetName("ProductionCuts");
387 AliRDHFCutsDplustoKpipi *analysis = new AliRDHFCutsDplustoKpipi(*fRDCutsAnalysis);
388 analysis->SetName("AnalysisCuts");
4c230f6d 389
390 fListCuts->Add(production);
391 fListCuts->Add(analysis);
392 PostData(2,fListCuts);
393
d486095a 394 return;
395}
396
397//________________________________________________________________________
398void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
399{
400 // Create the output container
401 //
402 if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");
403
404 // Several histograms are more conveniently managed in a TList
405 fOutput = new TList();
406 fOutput->SetOwner();
1f4e9722 407 fOutput->SetName("OutputHistos");
408
1f4e9722 409 TString hisname;
82bb8d4b 410 Int_t index=0;
993bfba1 411 Int_t nbins=GetNBinsHistos();
37fc80e6 412 fHistCentrality[0]=new TH1F("centrality","centrality",100,0.5,30000.5);
413 fHistCentrality[1]=new TH1F("centrality(selectedCent)","centrality(selectedCent)",100,0.5,30000.5);
414 fHistCentrality[2]=new TH1F("centrality(OutofCent)","centrality(OutofCent)",100,0.5,30000.5);
415 for(Int_t i=0;i<3;i++){
416 fHistCentrality[i]->Sumw2();
417 fOutput->Add(fHistCentrality[i]);
418 }
82bb8d4b 419 for(Int_t i=0;i<fNPtBins;i++){
1f4e9722 420
82bb8d4b 421 index=GetHistoIndex(i);
ba9ae5b2 422
82bb8d4b 423 hisname.Form("hMassPt%d",i);
993bfba1 424 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 425 fMassHist[index]->Sumw2();
ba9ae5b2 426 hisname.Form("hCosPAllPt%d",i);
993bfba1 427 fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
ba9ae5b2 428 fCosPHist[index]->Sumw2();
429 hisname.Form("hDLenAllPt%d",i);
993bfba1 430 fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
ba9ae5b2 431 fDLenHist[index]->Sumw2();
432 hisname.Form("hSumd02AllPt%d",i);
993bfba1 433 fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
ba9ae5b2 434 fSumd02Hist[index]->Sumw2();
435 hisname.Form("hSigVertAllPt%d",i);
993bfba1 436 fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
ba9ae5b2 437 fSigVertHist[index]->Sumw2();
438 hisname.Form("hPtMaxAllPt%d",i);
993bfba1 439 fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
ba9ae5b2 440 fPtMaxHist[index]->Sumw2();
3ec9254b 441 hisname.Form("hPtKPt%d",i);
442 fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
443 fPtKHist[index]->Sumw2();
444 hisname.Form("hPtpi1Pt%d",i);
445 fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
446 fPtpi1Hist[index]->Sumw2();
447 hisname.Form("hPtpi2Pt%d",i);
448 fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
449 fPtpi2Hist[index]->Sumw2();
ba9ae5b2 450 hisname.Form("hDCAAllPt%d",i);
993bfba1 451 fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
ba9ae5b2 452 fDCAHist[index]->Sumw2();
453
3ec9254b 454 hisname.Form("hDLxyPt%d",i);
455 fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
456 fDLxy[index]->Sumw2();
457 hisname.Form("hCosxyPt%d",i);
458 fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
459 fCosxy[index]->Sumw2();
460 hisname.Form("hDLxyPt%dTC",i);
461 fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
462 fDLxyTC[index]->Sumw2();
463 hisname.Form("hCosxyPt%dTC",i);
464 fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
465 fCosxyTC[index]->Sumw2();
466
1f4e9722 467 hisname.Form("hMassPt%dTC",i);
993bfba1 468 fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 469 fMassHistTC[index]->Sumw2();
bb69f127 470 hisname.Form("hMassPt%dTCPlus",i);
471 fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
472 fMassHistTCPlus[index]->Sumw2();
473 hisname.Form("hMassPt%dTCMinus",i);
474 fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
475 fMassHistTCMinus[index]->Sumw2();
ba9ae5b2 476
477
7d9ca2b5 478
82bb8d4b 479 index=GetSignalHistoIndex(i);
82bb8d4b 480 hisname.Form("hSigPt%d",i);
993bfba1 481 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 482 fMassHist[index]->Sumw2();
ba9ae5b2 483 hisname.Form("hCosPSigPt%d",i);
993bfba1 484 fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
ba9ae5b2 485 fCosPHist[index]->Sumw2();
486 hisname.Form("hDLenSigPt%d",i);
993bfba1 487 fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
ba9ae5b2 488 fDLenHist[index]->Sumw2();
489 hisname.Form("hSumd02SigPt%d",i);
993bfba1 490 fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
ba9ae5b2 491 fSumd02Hist[index]->Sumw2();
492 hisname.Form("hSigVertSigPt%d",i);
993bfba1 493 fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
ba9ae5b2 494 fSigVertHist[index]->Sumw2();
495 hisname.Form("hPtMaxSigPt%d",i);
993bfba1 496 fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
3ec9254b 497 fPtMaxHist[index]->Sumw2();
498 hisname.Form("hPtKSigPt%d",i);
499 fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
500 fPtKHist[index]->Sumw2();
501 hisname.Form("hPtpi1SigPt%d",i);
502 fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
503 fPtpi1Hist[index]->Sumw2();
504 hisname.Form("hPtpi2SigPt%d",i);
505 fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
506 fPtpi2Hist[index]->Sumw2();
ba9ae5b2 507
508 hisname.Form("hDCASigPt%d",i);
993bfba1 509 fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
ba9ae5b2 510 fDCAHist[index]->Sumw2();
511
3ec9254b 512 hisname.Form("hDLxySigPt%d",i);
513 fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
514 fDLxy[index]->Sumw2();
515 hisname.Form("hCosxySigPt%d",i);
516 fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
517 fCosxy[index]->Sumw2();
518 hisname.Form("hDLxySigPt%dTC",i);
519 fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
520 fDLxyTC[index]->Sumw2();
521 hisname.Form("hCosxySigPt%dTC",i);
522 fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
523 fCosxyTC[index]->Sumw2();
1f4e9722 524 hisname.Form("hSigPt%dTC",i);
993bfba1 525 fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 526 fMassHistTC[index]->Sumw2();
bb69f127 527 hisname.Form("hSigPt%dTCPlus",i);
528 fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
529 fMassHistTCPlus[index]->Sumw2();
530 hisname.Form("hSigPt%dTCMinus",i);
531 fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
532 fMassHistTCMinus[index]->Sumw2();
ba9ae5b2 533
ba9ae5b2 534
82bb8d4b 535 index=GetBackgroundHistoIndex(i);
82bb8d4b 536 hisname.Form("hBkgPt%d",i);
993bfba1 537 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 538 fMassHist[index]->Sumw2();
ba9ae5b2 539 hisname.Form("hCosPBkgPt%d",i);
993bfba1 540 fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
ba9ae5b2 541 fCosPHist[index]->Sumw2();
542 hisname.Form("hDLenBkgPt%d",i);
993bfba1 543 fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
ba9ae5b2 544 fDLenHist[index]->Sumw2();
545 hisname.Form("hSumd02BkgPt%d",i);
993bfba1 546 fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
ba9ae5b2 547 fSumd02Hist[index]->Sumw2();
548 hisname.Form("hSigVertBkgPt%d",i);
993bfba1 549 fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
ba9ae5b2 550 fSigVertHist[index]->Sumw2();
551 hisname.Form("hPtMaxBkgPt%d",i);
993bfba1 552 fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
ba9ae5b2 553 fPtMaxHist[index]->Sumw2();
3ec9254b 554 hisname.Form("hPtKBkgPt%d",i);
555 fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
556 fPtKHist[index]->Sumw2();
557 hisname.Form("hPtpi1BkgPt%d",i);
558 fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
559 fPtpi1Hist[index]->Sumw2();
560 hisname.Form("hPtpi2BkgPt%d",i);
561 fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
562 fPtpi2Hist[index]->Sumw2();
ba9ae5b2 563 hisname.Form("hDCABkgPt%d",i);
993bfba1 564 fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
ba9ae5b2 565 fDCAHist[index]->Sumw2();
566
3ec9254b 567 hisname.Form("hDLxyBkgPt%d",i);
568 fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
569 fDLxy[index]->Sumw2();
570 hisname.Form("hCosxyBkgPt%d",i);
571 fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
572 fCosxy[index]->Sumw2();
573 hisname.Form("hDLxyBkgPt%dTC",i);
574 fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
575 fDLxyTC[index]->Sumw2();
576 hisname.Form("hCosxyBkgPt%dTC",i);
577 fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
578 fCosxyTC[index]->Sumw2();
579
ba9ae5b2 580
1f4e9722 581 hisname.Form("hBkgPt%dTC",i);
993bfba1 582 fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 583 fMassHistTC[index]->Sumw2();
bb69f127 584 hisname.Form("hBkgPt%dTCPlus",i);
585 fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
586 fMassHistTCPlus[index]->Sumw2();
587 hisname.Form("hBkgPt%dTCMinus",i);
588 fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
589 fMassHistTCMinus[index]->Sumw2();
1f4e9722 590 }
7d9ca2b5 591
3ec9254b 592
82bb8d4b 593 for(Int_t i=0; i<3*fNPtBins; i++){
594 fOutput->Add(fMassHist[i]);
2edd194d 595 if(fCutsDistr){
596 fOutput->Add(fCosPHist[i]);
597 fOutput->Add(fDLenHist[i]);
598 fOutput->Add(fSumd02Hist[i]);
599 fOutput->Add(fSigVertHist[i]);
600 fOutput->Add(fPtMaxHist[i]);
601 fOutput->Add(fPtKHist[i]);
602 fOutput->Add(fPtpi1Hist[i]);
603 fOutput->Add(fPtpi2Hist[i]);
604 fOutput->Add(fDCAHist[i]);
605 fOutput->Add(fDLxy[i]);
606 fOutput->Add(fDLxyTC[i]);
607 fOutput->Add(fCosxy[i]);
608 fOutput->Add(fCosxyTC[i]);
609 }
82bb8d4b 610 fOutput->Add(fMassHistTC[i]);
bb69f127 611 fOutput->Add(fMassHistTCPlus[i]);
612 fOutput->Add(fMassHistTCMinus[i]);
2edd194d 613
82bb8d4b 614 }
ba9ae5b2 615
cdc197b1 616 if(fCutsDistr){
617 fCorreld0Kd0pi[0]=new TH2F("hCorreld0Kd0piAll","",100,-0.02,0.02,100,-0.02,0.02);
618 fCorreld0Kd0pi[1]=new TH2F("hCorreld0Kd0piSig","",100,-0.02,0.02,100,-0.02,0.02);
619 fCorreld0Kd0pi[2]=new TH2F("hCorreld0Kd0piBkg","",100,-0.02,0.02,100,-0.02,0.02);
620 for(Int_t i=0; i<3; i++){
621 fCorreld0Kd0pi[i]->Sumw2();
622 fOutput->Add(fCorreld0Kd0pi[i]);
623 }
624 }
26da4c42 625
5fc4893f 626 fHistNEvents = new TH1F("fHistNEvents", "number of events ",9,-0.5,8.5);
de1b82ff 627 fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
26da4c42 628 fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with good vertex");
de1b82ff 629 fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents with PbPb HM trigger");
630 fHistNEvents->GetXaxis()->SetBinLabel(4,"no. of Rejected pileup events");
631 fHistNEvents->GetXaxis()->SetBinLabel(5,"no. of candidate");
632 fHistNEvents->GetXaxis()->SetBinLabel(6,"no. of D+ after loose cuts");
633 fHistNEvents->GetXaxis()->SetBinLabel(7,"no. of D+ after tight cuts");
37fc80e6 634 fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events");
5fc4893f 635 fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
7d9ca2b5 636 fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
ba9ae5b2 637 fHistNEvents->Sumw2();
638 fHistNEvents->SetMinimum(0);
639 fOutput->Add(fHistNEvents);
595cc7e2 640
4d797793 641 fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
642 fPtVsMassTC=new TH2F("hPtVsMassTC","PtVsMass (analysis cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
595cc7e2 643 fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
644 fYVsPtTC=new TH2F("hYVsPtTC","YvsPt (analysis cuts)",40,0.,20.,80,-2.,2.);
645 fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
646 fYVsPtSigTC=new TH2F("hYVsPtSigTC","YvsPt (MC, only Sig, analysis cuts)",40,0.,20.,80,-2.,2.);
647
648 fOutput->Add(fPtVsMass);
649 fOutput->Add(fPtVsMassTC);
650 fOutput->Add(fYVsPt);
651 fOutput->Add(fYVsPtTC);
652 fOutput->Add(fYVsPtSig);
653 fOutput->Add(fYVsPtSigTC);
4afc48a2 654
a96083b9 655
656 //Counter for Normalization
e0fdfa52 657 TString normName="NormalizationCounter";
658 AliAnalysisDataContainer *cont = GetOutputSlot(3)->GetContainer();
659 if(cont)normName=(TString)cont->GetName();
660 fCounter = new AliNormalizationCounter(normName.Data());
38802708 661 fCounter->Init();
a96083b9 662
7d9ca2b5 663 if(fDoLS) CreateLikeSignHistos();
664 if(fDoImpPar) CreateImpactParameterHistos();
665
1f4e9722 666 if(fFillNtuple){
e0fdfa52 667 OpenFile(4); // 4 is the slot number of the ntuple
ba9ae5b2 668
669 fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:PtTrue:VxTrue:VyTrue:VzTrue:Ptpi:PtK:Ptpi2:PtRec:PointingAngle:DecLeng:VxRec:VyRec:VzRec:InvMass:sigvert:d0Pi:d0K:d0Pi2:dca:d0square");
670
1f4e9722 671 }
ba9ae5b2 672
d486095a 673 return;
674}
675
676//________________________________________________________________________
677void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
678{
679 // Execute analysis for current event:
680 // heavy flavor candidates association to MC truth
681
4c230f6d 682 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
4afc48a2 683
a96083b9 684
685
b557eb43 686 TClonesArray *array3Prong = 0;
4afc48a2 687 TClonesArray *arrayLikeSign =0;
b557eb43 688 if(!aod && AODEvent() && IsStandardAOD()) {
689 // In case there is an AOD handler writing a standard AOD, use the AOD
690 // event in memory rather than the input (ESD) event.
691 aod = dynamic_cast<AliAODEvent*> (AODEvent());
4c230f6d 692 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
693 // have to taken from the AOD event hold by the AliAODExtension
b557eb43 694 AliAODHandler* aodHandler = (AliAODHandler*)
695 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
696 if(aodHandler->GetExtensions()) {
697 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
698 AliAODEvent *aodFromExt = ext->GetAOD();
699 array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
4afc48a2 700 arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
b557eb43 701 }
dc222f77 702 } else if(aod) {
b557eb43 703 array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
4afc48a2 704 arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
b557eb43 705 }
8931c313 706
90993728 707 if(!aod || !array3Prong) {
d486095a 708 printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
709 return;
710 }
3a3c90f4 711 if(!arrayLikeSign && fDoLS) {
4afc48a2 712 printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
82bb8d4b 713 return;
4afc48a2 714 }
715
7c23877d 716
717 // fix for temporary bug in ESDfilter
718 // the AODs with null vertex pointer didn't pass the PhysSel
a96083b9 719 if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
1879baec 720 fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC);
de1b82ff 721 fHistNEvents->Fill(0); // count event
37fc80e6 722
723 Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
724 Float_t centrality=aod->GetNTracks();//fRDCutsAnalysis->GetCentrality(aod);
725 fHistCentrality[0]->Fill(centrality);
de1b82ff 726 // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
727 TString trigclass=aod->GetFiredTriggerClasses();
728 if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")||trigclass.Contains("C0SMH-B-NOPF-ALL")) fHistNEvents->Fill(2);
729 if(fRDCutsAnalysis->GetWhyRejection()==1)fHistNEvents->Fill(3);
37fc80e6 730 if(fRDCutsAnalysis->GetWhyRejection()==2){fHistNEvents->Fill(7);fHistCentrality[2]->Fill(centrality);}
731
732 // Post the data already here
de1b82ff 733 PostData(1,fOutput);
37fc80e6 734 if(!isEvSel)return;
735
736 fHistCentrality[1]->Fill(centrality);
737 fHistNEvents->Fill(1);
738
4afc48a2 739 TClonesArray *arrayMC=0;
740 AliAODMCHeader *mcHeader=0;
d486095a 741
742 // AOD primary vertex
1f4e9722 743 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
744 // vtx1->Print();
26da4c42 745 TString primTitle = vtx1->GetTitle();
746 //if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0)fHistNEvents->Fill(2);
747
d486095a 748 // load MC particles
4afc48a2 749 if(fReadMC){
750
751 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
752 if(!arrayMC) {
753 printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
90993728 754 return;
4afc48a2 755 }
1f4e9722 756
d486095a 757 // load MC header
4afc48a2 758 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
759 if(!mcHeader) {
82bb8d4b 760 printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
761 return;
4afc48a2 762 }
d486095a 763 }
4afc48a2 764
1f4e9722 765 Int_t n3Prong = array3Prong->GetEntriesFast();
7d9ca2b5 766 // printf("Number of D+->Kpipi: %d and of tracks: %d\n",n3Prong,aod->GetNumberOfTracks());
4afc48a2 767
768
769 Int_t nOS=0;
82bb8d4b 770 Int_t index;
1f4e9722 771 Int_t pdgDgDplustoKpipi[3]={321,211,211};
3ec9254b 772
773 if(fDoLS>1){
774 for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
775 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
776 if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
777 if(fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod))nOS++;
778 }
779 }
780 }else{
4c230f6d 781 // 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
782 //Double_t *cutsDplus = new (Double_t*)fRDCuts->GetCuts();
a96083b9 783 Int_t nSelectedloose=0,nSelectedtight=0;
1f4e9722 784 for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
785 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
de1b82ff 786 fHistNEvents->Fill(4);
5fc4893f 787 if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
788 fHistNEvents->Fill(8);
789 continue;
790 }
1f4e9722 791 Bool_t unsetvtx=kFALSE;
792 if(!d->GetOwnPrimaryVtx()){
793 d->SetOwnPrimaryVtx(vtx1);
794 unsetvtx=kTRUE;
795 }
4afc48a2 796
1275bf81 797 if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)) {
a96083b9 798
1f4e9722 799 Double_t ptCand = d->Pt();
a6003e0a 800 Int_t iPtBin = fRDCutsProduction->PtBin(ptCand);
82bb8d4b 801
1275bf81 802 Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
a6003e0a 803 Bool_t recVtx=kFALSE;
804 AliAODVertex *origownvtx=0x0;
805 if(fRDCutsProduction->GetIsPrimaryWithoutDaughters()){
806 if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
807 if(fRDCutsProduction->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
808 else fRDCutsProduction->CleanOwnPrimaryVtx(d,aod,origownvtx);
809 }
2bf2e62b 810
4afc48a2 811 Int_t labDp=-1;
7d9ca2b5 812 Bool_t isPrimary=kTRUE;
4afc48a2 813 Float_t deltaPx=0.;
814 Float_t deltaPy=0.;
815 Float_t deltaPz=0.;
816 Float_t truePt=0.;
817 Float_t xDecay=0.;
818 Float_t yDecay=0.;
819 Float_t zDecay=0.;
820 Float_t pdgCode=-2;
af636507 821 Float_t trueImpParXY=0.;
4afc48a2 822 if(fReadMC){
823 labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
824 if(labDp>=0){
825 AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
7d9ca2b5 826 if(CheckOrigin(arrayMC,partDp)==5) isPrimary=kFALSE;
4afc48a2 827 AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
828 deltaPx=partDp->Px()-d->Px();
829 deltaPy=partDp->Py()-d->Py();
830 deltaPz=partDp->Pz()-d->Pz();
831 truePt=partDp->Pt();
832 xDecay=dg0->Xv();
833 yDecay=dg0->Yv();
834 zDecay=dg0->Zv();
835 pdgCode=TMath::Abs(partDp->GetPdgCode());
af636507 836 if(!isPrimary){
94da94c4 837 trueImpParXY=GetTrueImpactParameter(mcHeader,arrayMC,partDp)*10000.;
af636507 838 }
4afc48a2 839 }else{
840 pdgCode=-1;
841 }
842 }
a6003e0a 843
1f4e9722 844 Double_t invMass=d->InvMassDplus();
595cc7e2 845 Double_t rapid=d->YDplus();
846 fYVsPt->Fill(ptCand,rapid);
3ec9254b 847 if(passTightCuts) {fYVsPtTC->Fill(ptCand,rapid);nOS++;}
595cc7e2 848 Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
849 if(isFidAcc){
850 fPtVsMass->Fill(invMass,ptCand);
851 if(passTightCuts) fPtVsMassTC->Fill(invMass,ptCand);
852 }
ba9ae5b2 853 Float_t tmp[24];
a6003e0a 854 Double_t dlen=d->DecayLength();
855 Double_t cosp=d->CosPointingAngle();
856 Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
857 Double_t dca=d->GetDCA();
858 Double_t sigvert=d->GetSigmaVert();
859 Double_t ptmax=0;
860 for(Int_t i=0;i<3;i++){
861 if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
862 }
7d9ca2b5 863 Double_t impparXY=d->ImpParXY()*10000.;
3de19caa 864 Double_t arrayForSparse[3]={invMass,ptCand,impparXY};
865 Double_t arrayForSparseTrue[3]={invMass,ptCand,trueImpParXY};
7d9ca2b5 866 if(fFillNtuple){
4afc48a2 867 tmp[0]=pdgCode;
868 tmp[1]=deltaPx;
869 tmp[2]=deltaPy;
870 tmp[3]=deltaPz;
871 tmp[4]=truePt;
872 tmp[5]=xDecay;
873 tmp[6]=yDecay;
874 tmp[7]=zDecay;
875 tmp[8]=d->PtProng(0);
876 tmp[9]=d->PtProng(1);
877 tmp[10]=d->PtProng(2);
878 tmp[11]=d->Pt();
a6003e0a 879 tmp[12]=cosp;
880 tmp[13]=dlen;
4afc48a2 881 tmp[14]=d->Xv();
882 tmp[15]=d->Yv();
883 tmp[16]=d->Zv();
884 tmp[17]=d->InvMassDplus();
a6003e0a 885 tmp[18]=sigvert;
4afc48a2 886 tmp[19]=d->Getd0Prong(0);
887 tmp[20]=d->Getd0Prong(1);
ba9ae5b2 888 tmp[21]=d->Getd0Prong(2);
a6003e0a 889 tmp[22]=dca;
ba9ae5b2 890 tmp[23]=d->Prodd0d0();
4afc48a2 891 fNtupleDplus->Fill(tmp);
de1b82ff 892 PostData(4,fNtupleDplus);
4afc48a2 893 }
82bb8d4b 894 if(iPtBin>=0){
3ec9254b 895 Float_t dlxy=d->NormalizedDecayLengthXY();
896 Float_t cxy=d->CosPointingAngleXY();
82bb8d4b 897 index=GetHistoIndex(iPtBin);
595cc7e2 898 if(isFidAcc){
de1b82ff 899 fHistNEvents->Fill(5);
a96083b9 900 nSelectedloose++;
595cc7e2 901 fMassHist[index]->Fill(invMass);
2edd194d 902 if(fCutsDistr){
903 fCosPHist[index]->Fill(cosp);
904 fDLenHist[index]->Fill(dlen);
905 fSumd02Hist[index]->Fill(sumD02);
906 fSigVertHist[index]->Fill(sigvert);
907 fPtMaxHist[index]->Fill(ptmax);
908 fPtKHist[index]->Fill(d->PtProng(1));
909 fPtpi1Hist[index]->Fill(d->PtProng(0));
910 fPtpi2Hist[index]->Fill(d->PtProng(2));
911 fDCAHist[index]->Fill(dca);
912 fDLxy[index]->Fill(dlxy);
913 fCosxy[index]->Fill(cxy);
cdc197b1 914 fCorreld0Kd0pi[0]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
915 d->Getd0Prong(2)*d->Getd0Prong(1));
2edd194d 916 }
de1b82ff 917 if(passTightCuts){ fHistNEvents->Fill(6);
a96083b9 918 nSelectedtight++;
595cc7e2 919 fMassHistTC[index]->Fill(invMass);
2edd194d 920 if(fCutsDistr){
921 fDLxyTC[index]->Fill(dlxy);
922 fCosxyTC[index]->Fill(cxy);
923 }
bb69f127 924 if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
925 else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
7d9ca2b5 926 if(fDoImpPar){
3de19caa 927 fHistMassPtImpParTC[0]->Fill(arrayForSparse);
7d9ca2b5 928 }
929 }
82bb8d4b 930 }
3ec9254b 931
82bb8d4b 932 if(fReadMC){
2edd194d 933 // if(fCutsDistr){
82bb8d4b 934 if(labDp>=0) {
935 index=GetSignalHistoIndex(iPtBin);
595cc7e2 936 if(isFidAcc){
a96083b9 937 Float_t factor[3]={1.,1.,1.};
938 if(fUseStrangeness){
939 for(Int_t iprong=0;iprong<3;iprong++){
940 AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
941 Int_t labd= trad->GetLabel();
942 if(labd>=0){
943 AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
944 if(dau){
945 Int_t labm = dau->GetMother();
946 if(labm>=0){
947 AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
948 if(mot){
949 if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
950 if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
951 else factor[iprong]=1./.6;
952 // fNentries->Fill(12);
953 }
954 if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
955 factor[iprong]=1./0.25;
956 // fNentries->Fill(13);
957 }//if 3122
958 }//if(mot)
959 }//if labm>0
960 }//if(dau)
961 }//if labd>=0
962 }//prong loop
963 }
964 Float_t fact=1.;for(Int_t k=0;k<3;k++)fact=fact*factor[k];
595cc7e2 965 fMassHist[index]->Fill(invMass);
2edd194d 966 if(fCutsDistr){
967 fCosPHist[index]->Fill(cosp,fact);
968 fDLenHist[index]->Fill(dlen,fact);
969 fDLxy[index]->Fill(dlxy);
970 fCosxy[index]->Fill(cxy);
3ec9254b 971
2edd194d 972 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];
973 fSumd02Hist[index]->Fill(sumd02s);
974 fSigVertHist[index]->Fill(sigvert,fact);
975 fPtMaxHist[index]->Fill(ptmax,fact);
976 fPtKHist[index]->Fill(d->PtProng(1),fact);
977 fPtpi1Hist[index]->Fill(d->PtProng(0),fact);
978 fPtpi2Hist[index]->Fill(d->PtProng(2),fact);
979 fDCAHist[index]->Fill(dca,fact);
cdc197b1 980 fCorreld0Kd0pi[1]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
981 d->Getd0Prong(2)*d->Getd0Prong(1));
2edd194d 982 }
595cc7e2 983 if(passTightCuts){
984 fMassHistTC[index]->Fill(invMass);
2edd194d 985 if(fCutsDistr){
986 fDLxyTC[index]->Fill(dlxy);
987 fCosxyTC[index]->Fill(cxy);
988 }
bb69f127 989 if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
990 else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
af636507 991 if(fDoImpPar){
3de19caa 992 if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse);
af636507 993 else{
3de19caa 994 fHistMassPtImpParTC[2]->Fill(arrayForSparse);
995 fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue);
af636507 996 }
997 }
595cc7e2 998 }
999 }
1000 fYVsPtSig->Fill(ptCand,rapid);
1001 if(passTightCuts) fYVsPtSigTC->Fill(ptCand,rapid);
82bb8d4b 1002 }else{
1003 index=GetBackgroundHistoIndex(iPtBin);
595cc7e2 1004 if(isFidAcc){
a96083b9 1005 Float_t factor[3]={1.,1.,1.};
1006 if(fUseStrangeness){
1007 for(Int_t iprong=0;iprong<3;iprong++){
1008 AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
1009 Int_t labd= trad->GetLabel();
1010 if(labd>=0){
1011 AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
1012 if(dau){
1013 Int_t labm = dau->GetMother();
1014 if(labm>=0){
1015 AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
1016 if(mot){
1017 if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
1018 if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
1019 else factor[iprong]=1./.6;
1020 // fNentries->Fill(12);
1021 }
1022 if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
1023 factor[iprong]=1./0.25;
1024 // fNentries->Fill(13);
1025 }//if 3122
1026 }//if(mot)
1027 }//if labm>0
1028 }//if(dau)
1029 }//if labd>=0
1030 }//prong loop
1031 }
2edd194d 1032
a96083b9 1033 Float_t fact=1.;for(Int_t k=0;k<3;k++)fact=fact*factor[k];
595cc7e2 1034 fMassHist[index]->Fill(invMass);
2edd194d 1035 if(fCutsDistr){
1036 fCosPHist[index]->Fill(cosp,fact);
1037 fDLenHist[index]->Fill(dlen,fact);
1038 fDLxy[index]->Fill(dlxy);
1039 fCosxy[index]->Fill(cxy);
3ec9254b 1040
2edd194d 1041 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];
1042 fSumd02Hist[index]->Fill(sumd02s);
1043 fSigVertHist[index]->Fill(sigvert,fact);
1044 fPtMaxHist[index]->Fill(ptmax,fact);
1045 fPtKHist[index]->Fill(d->PtProng(1),fact);
1046 fPtpi1Hist[index]->Fill(d->PtProng(0),fact);
1047 fPtpi2Hist[index]->Fill(d->PtProng(2),fact);
1048 fDCAHist[index]->Fill(dca,fact);
cdc197b1 1049 fCorreld0Kd0pi[2]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
1050 d->Getd0Prong(2)*d->Getd0Prong(1));
2edd194d 1051 }
595cc7e2 1052 if(passTightCuts){
1053 fMassHistTC[index]->Fill(invMass);
2edd194d 1054 if(fCutsDistr){
1055 fDLxyTC[index]->Fill(dlxy);
1056 fCosxyTC[index]->Fill(cxy);
1057 }
bb69f127 1058 if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
1059 else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
af636507 1060 if(fDoImpPar){
3de19caa 1061 fHistMassPtImpParTC[4]->Fill(arrayForSparse);
af636507 1062 }
595cc7e2 1063 }
3ec9254b 1064 }
2bf2e62b 1065 }
3ec9254b 1066
1067 }
3ec9254b 1068 }
2edd194d 1069
a6003e0a 1070 if(recVtx)fRDCutsProduction->CleanOwnPrimaryVtx(d,aod,origownvtx);
d486095a 1071 }
1f4e9722 1072 if(unsetvtx) d->UnsetOwnPrimaryVtx();
1073 }
a96083b9 1074 fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
1075 fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
3ec9254b 1076 }
4afc48a2 1077 //start LS analysis
1078 if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
82bb8d4b 1079
ad42e35b 1080 PostData(1,fOutput);
1081 PostData(2,fListCuts);
a96083b9 1082 PostData(3,fCounter);
d486095a 1083 return;
1084}
1085
7d9ca2b5 1086//________________________________________________________________________
1087void AliAnalysisTaskSEDplus::CreateLikeSignHistos(){
1088 // Histos for Like Sign bckground
1089
1090 TString hisname;
1091 Int_t indexLS=0;
1092 Int_t index=0;
1093 Int_t nbins=GetNBinsHistos();
1094 for(Int_t i=0;i<fNPtBins;i++){
1095
1096 index=GetHistoIndex(i);
1097 indexLS=GetLSHistoIndex(i);
1098
1099 hisname.Form("hLSPt%dLC",i);
1100 fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1101 fMassHistLS[indexLS]->Sumw2();
1102 hisname.Form("hLSPt%dTC",i);
1103 fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1104 fMassHistLSTC[indexLS]->Sumw2();
1105
1106 hisname.Form("hCosPAllPt%dLS",i);
1107 fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1108 fCosPHistLS[index]->Sumw2();
1109 hisname.Form("hDLenAllPt%dLS",i);
1110 fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1111 fDLenHistLS[index]->Sumw2();
1112 hisname.Form("hSumd02AllPt%dLS",i);
1113 fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1114 fSumd02HistLS[index]->Sumw2();
1115 hisname.Form("hSigVertAllPt%dLS",i);
1116 fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1117 fSigVertHistLS[index]->Sumw2();
1118 hisname.Form("hPtMaxAllPt%dLS",i);
1119 fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1120 fPtMaxHistLS[index]->Sumw2();
1121 hisname.Form("hDCAAllPt%dLS",i);
1122 fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1123 fDCAHistLS[index]->Sumw2();
1124
1125 index=GetSignalHistoIndex(i);
1126 indexLS++;
1127
1128 hisname.Form("hLSPt%dLCnw",i);
1129 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1130 fMassHistLS[indexLS]->Sumw2();
1131 hisname.Form("hLSPt%dTCnw",i);
1132 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1133 fMassHistLSTC[indexLS]->Sumw2();
1134
1135 hisname.Form("hCosPSigPt%dLS",i);
1136 fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1137 fCosPHistLS[index]->Sumw2();
1138 hisname.Form("hDLenSigPt%dLS",i);
1139 fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1140 fDLenHistLS[index]->Sumw2();
1141 hisname.Form("hSumd02SigPt%dLS",i);
1142 fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1143 fSumd02HistLS[index]->Sumw2();
1144 hisname.Form("hSigVertSigPt%dLS",i);
1145 fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1146 fSigVertHistLS[index]->Sumw2();
1147 hisname.Form("hPtMaxSigPt%dLS",i);
1148 fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1149 fPtMaxHistLS[index]->Sumw2();
1150 hisname.Form("hDCASigPt%dLS",i);
1151 fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1152 fDCAHistLS[index]->Sumw2();
1153
1154 index=GetBackgroundHistoIndex(i);
1155 indexLS++;
1156
1157 hisname.Form("hLSPt%dLCntrip",i);
1158 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1159 fMassHistLS[indexLS]->Sumw2();
1160 hisname.Form("hLSPt%dTCntrip",i);
1161 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1162 fMassHistLSTC[indexLS]->Sumw2();
1163
1164 hisname.Form("hCosPBkgPt%dLS",i);
1165 fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1166 fCosPHistLS[index]->Sumw2();
1167 hisname.Form("hDLenBkgPt%dLS",i);
1168 fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1169 fDLenHistLS[index]->Sumw2();
1170 hisname.Form("hSumd02BkgPt%dLS",i);
1171 fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1172 fSumd02HistLS[index]->Sumw2();
1173 hisname.Form("hSigVertBkgPt%dLS",i);
1174 fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1175 fSigVertHistLS[index]->Sumw2();
1176 hisname.Form("hPtMaxBkgPt%dLS",i);
1177 fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1178 fPtMaxHistLS[index]->Sumw2();
1179 hisname.Form("hDCABkgPt%dLS",i);
1180 fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1181 fDCAHistLS[index]->Sumw2();
1182
1183 indexLS++;
1184 hisname.Form("hLSPt%dLCntripsinglecut",i);
1185 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1186 fMassHistLS[indexLS]->Sumw2();
1187 hisname.Form("hLSPt%dTCntripsinglecut",i);
1188 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1189 fMassHistLSTC[indexLS]->Sumw2();
1190
1191 indexLS++;
1192 hisname.Form("hLSPt%dLCspc",i);
1193 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1194 fMassHistLS[indexLS]->Sumw2();
1195 hisname.Form("hLSPt%dTCspc",i);
1196 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1197 fMassHistLSTC[indexLS]->Sumw2();
1198 }
1199
1200 for(Int_t i=0; i<3*fNPtBins; i++){
1201 fOutput->Add(fCosPHistLS[i]);
1202 fOutput->Add(fDLenHistLS[i]);
1203 fOutput->Add(fSumd02HistLS[i]);
1204 fOutput->Add(fSigVertHistLS[i]);
1205 fOutput->Add(fPtMaxHistLS[i]);
1206 fOutput->Add(fDCAHistLS[i]);
1207 }
1208 for(Int_t i=0; i<5*fNPtBins; i++){
1209 fOutput->Add(fMassHistLS[i]);
1210 fOutput->Add(fMassHistLSTC[i]);
1211 }
1212}
1213
1214//________________________________________________________________________
1215void AliAnalysisTaskSEDplus::CreateImpactParameterHistos(){
1216 // Histos for impact paramter study
4afc48a2 1217
3de19caa 1218 Int_t nmassbins=GetNBinsHistos();
1219 Int_t nbins[3]={nmassbins,200,fNImpParBins};
1220 Double_t xmin[3]={fLowmasslimit,0.,fLowerImpPar};
1221 Double_t xmax[3]={fUpmasslimit,20.,fHigherImpPar};
1222
1223 fHistMassPtImpParTC[0]=new THnSparseF("hMassPtImpParAll",
1224 "Mass vs. pt vs.imppar - All",
1225 3,nbins,xmin,xmax);
1226 fHistMassPtImpParTC[1]=new THnSparseF("hMassPtImpParPrompt",
1227 "Mass vs. pt vs.imppar - promptD",
1228 3,nbins,xmin,xmax);
1229 fHistMassPtImpParTC[2]=new THnSparseF("hMassPtImpParBfeed",
1230 "Mass vs. pt vs.imppar - DfromB",
1231 3,nbins,xmin,xmax);
1232 fHistMassPtImpParTC[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1233 "Mass vs. pt vs.true imppar -DfromB",
1234 3,nbins,xmin,xmax);
1235 fHistMassPtImpParTC[4]=new THnSparseF("hMassPtImpParBkg",
1236 "Mass vs. pt vs.imppar - backgr.",
1237 3,nbins,xmin,xmax);
1238
7d9ca2b5 1239 for(Int_t i=0; i<5;i++){
1240 fOutput->Add(fHistMassPtImpParTC[i]);
1241 }
1242}
4afc48a2 1243
82bb8d4b 1244//________________________________________________________________________
d486095a 1245void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
1246{
1247 // Terminate analysis
1248 //
1249 if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
82bb8d4b 1250
d486095a 1251 fOutput = dynamic_cast<TList*> (GetOutputData(1));
1252 if (!fOutput) {
1253 printf("ERROR: fOutput not available\n");
1254 return;
1255 }
595cc7e2 1256 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
82bb8d4b 1257
595cc7e2 1258 TString hisname;
1259 Int_t index=0;
ba9ae5b2 1260
595cc7e2 1261 for(Int_t i=0;i<fNPtBins;i++){
82bb8d4b 1262 index=GetHistoIndex(i);
ba9ae5b2 1263
37fc80e6 1264 hisname.Form("hMassPt%dTC",i);
82bb8d4b 1265 fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
37fc80e6 1266 }
ba9ae5b2 1267
13808a30 1268 TCanvas *c1=new TCanvas("c1","D+ invariant mass distribution",500,500);
1269 c1->cd();
1270 TH1F *hMassPt=(TH1F*)fOutput->FindObject("hMassPt3TC");
1271 hMassPt->SetLineColor(kBlue);
1272 hMassPt->SetXTitle("M[GeV/c^{2}]");
1273 hMassPt->Draw();
1274
595cc7e2 1275 return;
d486095a 1276}
7d9ca2b5 1277//_________________________________________________________________________________________________
1278Int_t AliAnalysisTaskSEDplus::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
1279 //
1280 // checking whether the mother of the particles come from a charm or a bottom quark
1281 //
1282
1283 Int_t pdgGranma = 0;
1284 Int_t mother = 0;
1285 mother = mcPartCandidate->GetMother();
1286 Int_t istep = 0;
1287 Int_t abspdgGranma =0;
1288 Bool_t isFromB=kFALSE;
1289 Bool_t isQuarkFound=kFALSE;
1290 while (mother >0 ){
1291 istep++;
1292 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1293 if (mcGranma){
1294 pdgGranma = mcGranma->GetPdgCode();
1295 abspdgGranma = TMath::Abs(pdgGranma);
1296 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1297 isFromB=kTRUE;
1298 }
1299 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1300 mother = mcGranma->GetMother();
1301 }else{
1302 AliError("Failed casting the mother particle!");
1303 break;
1304 }
1305 }
1306
1307 if(isFromB) return 5;
1308 else return 4;
1309}
af636507 1310//_________________________________________________________________________________________________
1311Float_t AliAnalysisTaskSEDplus::GetTrueImpactParameter(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partDp) const {
1312 // true impact parameter calculation
1313
1314 Double_t vtxTrue[3];
1315 mcHeader->GetVertex(vtxTrue);
1316 Double_t origD[3];
94da94c4 1317 partDp->XvYvZv(origD);
af636507 1318 Short_t charge=partDp->Charge();
1319 Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
1320 for(Int_t iDau=0; iDau<3; iDau++){
1321 pXdauTrue[iDau]=0.;
1322 pYdauTrue[iDau]=0.;
1323 pZdauTrue[iDau]=0.;
1324 }
1325
1326 Int_t nDau=partDp->GetNDaughters();
1327 Int_t labelFirstDau = partDp->GetDaughter(0);
1328 if(nDau==3){
1329 for(Int_t iDau=0; iDau<3; iDau++){
1330 Int_t ind = labelFirstDau+iDau;
1331 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
3de19caa 1332 if(!part){
1333 AliError("Daughter particle not found in MC array");
1334 return 99999.;
1335 }
af636507 1336 pXdauTrue[iDau]=part->Px();
1337 pYdauTrue[iDau]=part->Py();
1338 pZdauTrue[iDau]=part->Pz();
1339 }
1340 }else if(nDau==2){
1341 Int_t theDau=0;
1342 for(Int_t iDau=0; iDau<2; iDau++){
1343 Int_t ind = labelFirstDau+iDau;
1344 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
3de19caa 1345 if(!part){
1346 AliError("Daughter particle not found in MC array");
1347 return 99999.;
1348 }
af636507 1349 Int_t pdgCode=TMath::Abs(part->GetPdgCode());
1350 if(pdgCode==211 || pdgCode==321){
1351 pXdauTrue[theDau]=part->Px();
1352 pYdauTrue[theDau]=part->Py();
1353 pZdauTrue[theDau]=part->Pz();
1354 ++theDau;
1355 }else{
1356 Int_t nDauRes=part->GetNDaughters();
1357 if(nDauRes==2){
1358 Int_t labelFirstDauRes = part->GetDaughter(0);
1359 for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
1360 Int_t indDR = labelFirstDauRes+iDauRes;
1361 AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
3de19caa 1362 if(!partDR){
1363 AliError("Daughter particle not found in MC array");
1364 return 99999.;
1365 }
1366
1367 Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
1368 if(pdgCodeDR==211 || pdgCodeDR==321){
af636507 1369 pXdauTrue[theDau]=partDR->Px();
1370 pYdauTrue[theDau]=partDR->Py();
1371 pZdauTrue[theDau]=partDR->Pz();
1372 ++theDau;
1373 }
1374 }
1375 }
1376 }
1377 }
3de19caa 1378 if(theDau!=3){
1379 AliError("Wrong number of decay prongs");
1380 return 99999.;
1381 }
af636507 1382 }
1383
1384 Double_t d0dummy[3]={0.,0.,0.};
1385 AliAODRecoDecayHF aodDplusMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
1386 return aodDplusMC.ImpParXY();
1387
1388}