]> 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
Coding conventions
[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);
c0d01d9d 724 Bool_t isEvSelP=kTRUE;
725 isEvSelP=fRDCutsProduction->IsEventSelected(aod); // to have proper PID object settings
97049c4a 726
37fc80e6 727 Float_t centrality=aod->GetNTracks();//fRDCutsAnalysis->GetCentrality(aod);
728 fHistCentrality[0]->Fill(centrality);
de1b82ff 729 // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
730 TString trigclass=aod->GetFiredTriggerClasses();
731 if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")||trigclass.Contains("C0SMH-B-NOPF-ALL")) fHistNEvents->Fill(2);
732 if(fRDCutsAnalysis->GetWhyRejection()==1)fHistNEvents->Fill(3);
37fc80e6 733 if(fRDCutsAnalysis->GetWhyRejection()==2){fHistNEvents->Fill(7);fHistCentrality[2]->Fill(centrality);}
734
735 // Post the data already here
de1b82ff 736 PostData(1,fOutput);
37fc80e6 737 if(!isEvSel)return;
738
739 fHistCentrality[1]->Fill(centrality);
740 fHistNEvents->Fill(1);
741
4afc48a2 742 TClonesArray *arrayMC=0;
743 AliAODMCHeader *mcHeader=0;
d486095a 744
745 // AOD primary vertex
1f4e9722 746 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
747 // vtx1->Print();
26da4c42 748 TString primTitle = vtx1->GetTitle();
749 //if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0)fHistNEvents->Fill(2);
750
d486095a 751 // load MC particles
4afc48a2 752 if(fReadMC){
753
754 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
755 if(!arrayMC) {
756 printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
90993728 757 return;
4afc48a2 758 }
1f4e9722 759
d486095a 760 // load MC header
4afc48a2 761 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
762 if(!mcHeader) {
82bb8d4b 763 printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
764 return;
4afc48a2 765 }
d486095a 766 }
4afc48a2 767
1f4e9722 768 Int_t n3Prong = array3Prong->GetEntriesFast();
7d9ca2b5 769 // printf("Number of D+->Kpipi: %d and of tracks: %d\n",n3Prong,aod->GetNumberOfTracks());
4afc48a2 770
771
772 Int_t nOS=0;
82bb8d4b 773 Int_t index;
1f4e9722 774 Int_t pdgDgDplustoKpipi[3]={321,211,211};
3ec9254b 775
776 if(fDoLS>1){
777 for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
778 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
779 if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
780 if(fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod))nOS++;
781 }
782 }
783 }else{
4c230f6d 784 // 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
785 //Double_t *cutsDplus = new (Double_t*)fRDCuts->GetCuts();
a96083b9 786 Int_t nSelectedloose=0,nSelectedtight=0;
1f4e9722 787 for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
788 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
de1b82ff 789 fHistNEvents->Fill(4);
5fc4893f 790 if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
791 fHistNEvents->Fill(8);
792 continue;
793 }
1f4e9722 794 Bool_t unsetvtx=kFALSE;
795 if(!d->GetOwnPrimaryVtx()){
796 d->SetOwnPrimaryVtx(vtx1);
797 unsetvtx=kTRUE;
798 }
4afc48a2 799
1275bf81 800 if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)) {
a96083b9 801
1f4e9722 802 Double_t ptCand = d->Pt();
a6003e0a 803 Int_t iPtBin = fRDCutsProduction->PtBin(ptCand);
82bb8d4b 804
1275bf81 805 Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
a6003e0a 806 Bool_t recVtx=kFALSE;
807 AliAODVertex *origownvtx=0x0;
808 if(fRDCutsProduction->GetIsPrimaryWithoutDaughters()){
809 if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
810 if(fRDCutsProduction->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
811 else fRDCutsProduction->CleanOwnPrimaryVtx(d,aod,origownvtx);
812 }
2bf2e62b 813
4afc48a2 814 Int_t labDp=-1;
7d9ca2b5 815 Bool_t isPrimary=kTRUE;
4afc48a2 816 Float_t deltaPx=0.;
817 Float_t deltaPy=0.;
818 Float_t deltaPz=0.;
819 Float_t truePt=0.;
820 Float_t xDecay=0.;
821 Float_t yDecay=0.;
822 Float_t zDecay=0.;
823 Float_t pdgCode=-2;
af636507 824 Float_t trueImpParXY=0.;
4afc48a2 825 if(fReadMC){
826 labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
827 if(labDp>=0){
828 AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
7d9ca2b5 829 if(CheckOrigin(arrayMC,partDp)==5) isPrimary=kFALSE;
4afc48a2 830 AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
831 deltaPx=partDp->Px()-d->Px();
832 deltaPy=partDp->Py()-d->Py();
833 deltaPz=partDp->Pz()-d->Pz();
834 truePt=partDp->Pt();
835 xDecay=dg0->Xv();
836 yDecay=dg0->Yv();
837 zDecay=dg0->Zv();
838 pdgCode=TMath::Abs(partDp->GetPdgCode());
af636507 839 if(!isPrimary){
94da94c4 840 trueImpParXY=GetTrueImpactParameter(mcHeader,arrayMC,partDp)*10000.;
af636507 841 }
4afc48a2 842 }else{
843 pdgCode=-1;
844 }
845 }
a6003e0a 846
1f4e9722 847 Double_t invMass=d->InvMassDplus();
595cc7e2 848 Double_t rapid=d->YDplus();
849 fYVsPt->Fill(ptCand,rapid);
3ec9254b 850 if(passTightCuts) {fYVsPtTC->Fill(ptCand,rapid);nOS++;}
595cc7e2 851 Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
852 if(isFidAcc){
853 fPtVsMass->Fill(invMass,ptCand);
854 if(passTightCuts) fPtVsMassTC->Fill(invMass,ptCand);
855 }
ba9ae5b2 856 Float_t tmp[24];
a6003e0a 857 Double_t dlen=d->DecayLength();
858 Double_t cosp=d->CosPointingAngle();
859 Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
860 Double_t dca=d->GetDCA();
861 Double_t sigvert=d->GetSigmaVert();
862 Double_t ptmax=0;
863 for(Int_t i=0;i<3;i++){
864 if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
865 }
7d9ca2b5 866 Double_t impparXY=d->ImpParXY()*10000.;
3de19caa 867 Double_t arrayForSparse[3]={invMass,ptCand,impparXY};
868 Double_t arrayForSparseTrue[3]={invMass,ptCand,trueImpParXY};
7d9ca2b5 869 if(fFillNtuple){
4afc48a2 870 tmp[0]=pdgCode;
871 tmp[1]=deltaPx;
872 tmp[2]=deltaPy;
873 tmp[3]=deltaPz;
874 tmp[4]=truePt;
875 tmp[5]=xDecay;
876 tmp[6]=yDecay;
877 tmp[7]=zDecay;
878 tmp[8]=d->PtProng(0);
879 tmp[9]=d->PtProng(1);
880 tmp[10]=d->PtProng(2);
881 tmp[11]=d->Pt();
a6003e0a 882 tmp[12]=cosp;
883 tmp[13]=dlen;
4afc48a2 884 tmp[14]=d->Xv();
885 tmp[15]=d->Yv();
886 tmp[16]=d->Zv();
887 tmp[17]=d->InvMassDplus();
a6003e0a 888 tmp[18]=sigvert;
4afc48a2 889 tmp[19]=d->Getd0Prong(0);
890 tmp[20]=d->Getd0Prong(1);
ba9ae5b2 891 tmp[21]=d->Getd0Prong(2);
a6003e0a 892 tmp[22]=dca;
ba9ae5b2 893 tmp[23]=d->Prodd0d0();
4afc48a2 894 fNtupleDplus->Fill(tmp);
de1b82ff 895 PostData(4,fNtupleDplus);
4afc48a2 896 }
82bb8d4b 897 if(iPtBin>=0){
3ec9254b 898 Float_t dlxy=d->NormalizedDecayLengthXY();
899 Float_t cxy=d->CosPointingAngleXY();
82bb8d4b 900 index=GetHistoIndex(iPtBin);
595cc7e2 901 if(isFidAcc){
de1b82ff 902 fHistNEvents->Fill(5);
a96083b9 903 nSelectedloose++;
595cc7e2 904 fMassHist[index]->Fill(invMass);
2edd194d 905 if(fCutsDistr){
906 fCosPHist[index]->Fill(cosp);
907 fDLenHist[index]->Fill(dlen);
908 fSumd02Hist[index]->Fill(sumD02);
909 fSigVertHist[index]->Fill(sigvert);
910 fPtMaxHist[index]->Fill(ptmax);
911 fPtKHist[index]->Fill(d->PtProng(1));
912 fPtpi1Hist[index]->Fill(d->PtProng(0));
913 fPtpi2Hist[index]->Fill(d->PtProng(2));
914 fDCAHist[index]->Fill(dca);
915 fDLxy[index]->Fill(dlxy);
916 fCosxy[index]->Fill(cxy);
cdc197b1 917 fCorreld0Kd0pi[0]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
918 d->Getd0Prong(2)*d->Getd0Prong(1));
2edd194d 919 }
de1b82ff 920 if(passTightCuts){ fHistNEvents->Fill(6);
a96083b9 921 nSelectedtight++;
595cc7e2 922 fMassHistTC[index]->Fill(invMass);
2edd194d 923 if(fCutsDistr){
924 fDLxyTC[index]->Fill(dlxy);
925 fCosxyTC[index]->Fill(cxy);
926 }
bb69f127 927 if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
928 else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
7d9ca2b5 929 if(fDoImpPar){
3de19caa 930 fHistMassPtImpParTC[0]->Fill(arrayForSparse);
7d9ca2b5 931 }
932 }
82bb8d4b 933 }
3ec9254b 934
82bb8d4b 935 if(fReadMC){
2edd194d 936 // if(fCutsDistr){
82bb8d4b 937 if(labDp>=0) {
938 index=GetSignalHistoIndex(iPtBin);
595cc7e2 939 if(isFidAcc){
a96083b9 940 Float_t factor[3]={1.,1.,1.};
941 if(fUseStrangeness){
942 for(Int_t iprong=0;iprong<3;iprong++){
943 AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
944 Int_t labd= trad->GetLabel();
945 if(labd>=0){
946 AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
947 if(dau){
948 Int_t labm = dau->GetMother();
949 if(labm>=0){
950 AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
951 if(mot){
952 if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
953 if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
954 else factor[iprong]=1./.6;
955 // fNentries->Fill(12);
956 }
957 if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
958 factor[iprong]=1./0.25;
959 // fNentries->Fill(13);
960 }//if 3122
961 }//if(mot)
962 }//if labm>0
963 }//if(dau)
964 }//if labd>=0
965 }//prong loop
966 }
967 Float_t fact=1.;for(Int_t k=0;k<3;k++)fact=fact*factor[k];
595cc7e2 968 fMassHist[index]->Fill(invMass);
2edd194d 969 if(fCutsDistr){
970 fCosPHist[index]->Fill(cosp,fact);
971 fDLenHist[index]->Fill(dlen,fact);
972 fDLxy[index]->Fill(dlxy);
973 fCosxy[index]->Fill(cxy);
3ec9254b 974
2edd194d 975 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];
976 fSumd02Hist[index]->Fill(sumd02s);
977 fSigVertHist[index]->Fill(sigvert,fact);
978 fPtMaxHist[index]->Fill(ptmax,fact);
979 fPtKHist[index]->Fill(d->PtProng(1),fact);
980 fPtpi1Hist[index]->Fill(d->PtProng(0),fact);
981 fPtpi2Hist[index]->Fill(d->PtProng(2),fact);
982 fDCAHist[index]->Fill(dca,fact);
cdc197b1 983 fCorreld0Kd0pi[1]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
984 d->Getd0Prong(2)*d->Getd0Prong(1));
2edd194d 985 }
595cc7e2 986 if(passTightCuts){
987 fMassHistTC[index]->Fill(invMass);
2edd194d 988 if(fCutsDistr){
989 fDLxyTC[index]->Fill(dlxy);
990 fCosxyTC[index]->Fill(cxy);
991 }
bb69f127 992 if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
993 else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
af636507 994 if(fDoImpPar){
3de19caa 995 if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse);
af636507 996 else{
3de19caa 997 fHistMassPtImpParTC[2]->Fill(arrayForSparse);
998 fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue);
af636507 999 }
1000 }
595cc7e2 1001 }
1002 }
1003 fYVsPtSig->Fill(ptCand,rapid);
1004 if(passTightCuts) fYVsPtSigTC->Fill(ptCand,rapid);
82bb8d4b 1005 }else{
1006 index=GetBackgroundHistoIndex(iPtBin);
595cc7e2 1007 if(isFidAcc){
a96083b9 1008 Float_t factor[3]={1.,1.,1.};
1009 if(fUseStrangeness){
1010 for(Int_t iprong=0;iprong<3;iprong++){
1011 AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
1012 Int_t labd= trad->GetLabel();
1013 if(labd>=0){
1014 AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
1015 if(dau){
1016 Int_t labm = dau->GetMother();
1017 if(labm>=0){
1018 AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
1019 if(mot){
1020 if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
1021 if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
1022 else factor[iprong]=1./.6;
1023 // fNentries->Fill(12);
1024 }
1025 if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
1026 factor[iprong]=1./0.25;
1027 // fNentries->Fill(13);
1028 }//if 3122
1029 }//if(mot)
1030 }//if labm>0
1031 }//if(dau)
1032 }//if labd>=0
1033 }//prong loop
1034 }
2edd194d 1035
a96083b9 1036 Float_t fact=1.;for(Int_t k=0;k<3;k++)fact=fact*factor[k];
595cc7e2 1037 fMassHist[index]->Fill(invMass);
2edd194d 1038 if(fCutsDistr){
1039 fCosPHist[index]->Fill(cosp,fact);
1040 fDLenHist[index]->Fill(dlen,fact);
1041 fDLxy[index]->Fill(dlxy);
1042 fCosxy[index]->Fill(cxy);
3ec9254b 1043
2edd194d 1044 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];
1045 fSumd02Hist[index]->Fill(sumd02s);
1046 fSigVertHist[index]->Fill(sigvert,fact);
1047 fPtMaxHist[index]->Fill(ptmax,fact);
1048 fPtKHist[index]->Fill(d->PtProng(1),fact);
1049 fPtpi1Hist[index]->Fill(d->PtProng(0),fact);
1050 fPtpi2Hist[index]->Fill(d->PtProng(2),fact);
1051 fDCAHist[index]->Fill(dca,fact);
cdc197b1 1052 fCorreld0Kd0pi[2]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
1053 d->Getd0Prong(2)*d->Getd0Prong(1));
2edd194d 1054 }
595cc7e2 1055 if(passTightCuts){
1056 fMassHistTC[index]->Fill(invMass);
2edd194d 1057 if(fCutsDistr){
1058 fDLxyTC[index]->Fill(dlxy);
1059 fCosxyTC[index]->Fill(cxy);
1060 }
bb69f127 1061 if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
1062 else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
af636507 1063 if(fDoImpPar){
3de19caa 1064 fHistMassPtImpParTC[4]->Fill(arrayForSparse);
af636507 1065 }
595cc7e2 1066 }
3ec9254b 1067 }
2bf2e62b 1068 }
3ec9254b 1069
1070 }
3ec9254b 1071 }
2edd194d 1072
a6003e0a 1073 if(recVtx)fRDCutsProduction->CleanOwnPrimaryVtx(d,aod,origownvtx);
d486095a 1074 }
1f4e9722 1075 if(unsetvtx) d->UnsetOwnPrimaryVtx();
1076 }
a96083b9 1077 fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
1078 fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
3ec9254b 1079 }
4afc48a2 1080 //start LS analysis
1081 if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
82bb8d4b 1082
ad42e35b 1083 PostData(1,fOutput);
1084 PostData(2,fListCuts);
a96083b9 1085 PostData(3,fCounter);
d486095a 1086 return;
1087}
1088
7d9ca2b5 1089//________________________________________________________________________
1090void AliAnalysisTaskSEDplus::CreateLikeSignHistos(){
1091 // Histos for Like Sign bckground
1092
1093 TString hisname;
1094 Int_t indexLS=0;
1095 Int_t index=0;
1096 Int_t nbins=GetNBinsHistos();
1097 for(Int_t i=0;i<fNPtBins;i++){
1098
1099 index=GetHistoIndex(i);
1100 indexLS=GetLSHistoIndex(i);
1101
1102 hisname.Form("hLSPt%dLC",i);
1103 fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1104 fMassHistLS[indexLS]->Sumw2();
1105 hisname.Form("hLSPt%dTC",i);
1106 fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1107 fMassHistLSTC[indexLS]->Sumw2();
1108
1109 hisname.Form("hCosPAllPt%dLS",i);
1110 fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1111 fCosPHistLS[index]->Sumw2();
1112 hisname.Form("hDLenAllPt%dLS",i);
1113 fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1114 fDLenHistLS[index]->Sumw2();
1115 hisname.Form("hSumd02AllPt%dLS",i);
1116 fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1117 fSumd02HistLS[index]->Sumw2();
1118 hisname.Form("hSigVertAllPt%dLS",i);
1119 fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1120 fSigVertHistLS[index]->Sumw2();
1121 hisname.Form("hPtMaxAllPt%dLS",i);
1122 fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1123 fPtMaxHistLS[index]->Sumw2();
1124 hisname.Form("hDCAAllPt%dLS",i);
1125 fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1126 fDCAHistLS[index]->Sumw2();
1127
1128 index=GetSignalHistoIndex(i);
1129 indexLS++;
1130
1131 hisname.Form("hLSPt%dLCnw",i);
1132 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1133 fMassHistLS[indexLS]->Sumw2();
1134 hisname.Form("hLSPt%dTCnw",i);
1135 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1136 fMassHistLSTC[indexLS]->Sumw2();
1137
1138 hisname.Form("hCosPSigPt%dLS",i);
1139 fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1140 fCosPHistLS[index]->Sumw2();
1141 hisname.Form("hDLenSigPt%dLS",i);
1142 fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1143 fDLenHistLS[index]->Sumw2();
1144 hisname.Form("hSumd02SigPt%dLS",i);
1145 fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1146 fSumd02HistLS[index]->Sumw2();
1147 hisname.Form("hSigVertSigPt%dLS",i);
1148 fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1149 fSigVertHistLS[index]->Sumw2();
1150 hisname.Form("hPtMaxSigPt%dLS",i);
1151 fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1152 fPtMaxHistLS[index]->Sumw2();
1153 hisname.Form("hDCASigPt%dLS",i);
1154 fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1155 fDCAHistLS[index]->Sumw2();
1156
1157 index=GetBackgroundHistoIndex(i);
1158 indexLS++;
1159
1160 hisname.Form("hLSPt%dLCntrip",i);
1161 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1162 fMassHistLS[indexLS]->Sumw2();
1163 hisname.Form("hLSPt%dTCntrip",i);
1164 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1165 fMassHistLSTC[indexLS]->Sumw2();
1166
1167 hisname.Form("hCosPBkgPt%dLS",i);
1168 fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1169 fCosPHistLS[index]->Sumw2();
1170 hisname.Form("hDLenBkgPt%dLS",i);
1171 fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1172 fDLenHistLS[index]->Sumw2();
1173 hisname.Form("hSumd02BkgPt%dLS",i);
1174 fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1175 fSumd02HistLS[index]->Sumw2();
1176 hisname.Form("hSigVertBkgPt%dLS",i);
1177 fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1178 fSigVertHistLS[index]->Sumw2();
1179 hisname.Form("hPtMaxBkgPt%dLS",i);
1180 fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1181 fPtMaxHistLS[index]->Sumw2();
1182 hisname.Form("hDCABkgPt%dLS",i);
1183 fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1184 fDCAHistLS[index]->Sumw2();
1185
1186 indexLS++;
1187 hisname.Form("hLSPt%dLCntripsinglecut",i);
1188 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1189 fMassHistLS[indexLS]->Sumw2();
1190 hisname.Form("hLSPt%dTCntripsinglecut",i);
1191 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1192 fMassHistLSTC[indexLS]->Sumw2();
1193
1194 indexLS++;
1195 hisname.Form("hLSPt%dLCspc",i);
1196 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1197 fMassHistLS[indexLS]->Sumw2();
1198 hisname.Form("hLSPt%dTCspc",i);
1199 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1200 fMassHistLSTC[indexLS]->Sumw2();
1201 }
1202
1203 for(Int_t i=0; i<3*fNPtBins; i++){
1204 fOutput->Add(fCosPHistLS[i]);
1205 fOutput->Add(fDLenHistLS[i]);
1206 fOutput->Add(fSumd02HistLS[i]);
1207 fOutput->Add(fSigVertHistLS[i]);
1208 fOutput->Add(fPtMaxHistLS[i]);
1209 fOutput->Add(fDCAHistLS[i]);
1210 }
1211 for(Int_t i=0; i<5*fNPtBins; i++){
1212 fOutput->Add(fMassHistLS[i]);
1213 fOutput->Add(fMassHistLSTC[i]);
1214 }
1215}
1216
1217//________________________________________________________________________
1218void AliAnalysisTaskSEDplus::CreateImpactParameterHistos(){
1219 // Histos for impact paramter study
4afc48a2 1220
3de19caa 1221 Int_t nmassbins=GetNBinsHistos();
1222 Int_t nbins[3]={nmassbins,200,fNImpParBins};
1223 Double_t xmin[3]={fLowmasslimit,0.,fLowerImpPar};
1224 Double_t xmax[3]={fUpmasslimit,20.,fHigherImpPar};
1225
1226 fHistMassPtImpParTC[0]=new THnSparseF("hMassPtImpParAll",
1227 "Mass vs. pt vs.imppar - All",
1228 3,nbins,xmin,xmax);
1229 fHistMassPtImpParTC[1]=new THnSparseF("hMassPtImpParPrompt",
1230 "Mass vs. pt vs.imppar - promptD",
1231 3,nbins,xmin,xmax);
1232 fHistMassPtImpParTC[2]=new THnSparseF("hMassPtImpParBfeed",
1233 "Mass vs. pt vs.imppar - DfromB",
1234 3,nbins,xmin,xmax);
1235 fHistMassPtImpParTC[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1236 "Mass vs. pt vs.true imppar -DfromB",
1237 3,nbins,xmin,xmax);
1238 fHistMassPtImpParTC[4]=new THnSparseF("hMassPtImpParBkg",
1239 "Mass vs. pt vs.imppar - backgr.",
1240 3,nbins,xmin,xmax);
1241
7d9ca2b5 1242 for(Int_t i=0; i<5;i++){
1243 fOutput->Add(fHistMassPtImpParTC[i]);
1244 }
1245}
4afc48a2 1246
82bb8d4b 1247//________________________________________________________________________
d486095a 1248void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
1249{
1250 // Terminate analysis
1251 //
1252 if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
82bb8d4b 1253
d486095a 1254 fOutput = dynamic_cast<TList*> (GetOutputData(1));
1255 if (!fOutput) {
1256 printf("ERROR: fOutput not available\n");
1257 return;
1258 }
595cc7e2 1259 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
82bb8d4b 1260
595cc7e2 1261 TString hisname;
1262 Int_t index=0;
ba9ae5b2 1263
595cc7e2 1264 for(Int_t i=0;i<fNPtBins;i++){
82bb8d4b 1265 index=GetHistoIndex(i);
ba9ae5b2 1266
37fc80e6 1267 hisname.Form("hMassPt%dTC",i);
82bb8d4b 1268 fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
37fc80e6 1269 }
ba9ae5b2 1270
13808a30 1271 TCanvas *c1=new TCanvas("c1","D+ invariant mass distribution",500,500);
1272 c1->cd();
1273 TH1F *hMassPt=(TH1F*)fOutput->FindObject("hMassPt3TC");
1274 hMassPt->SetLineColor(kBlue);
1275 hMassPt->SetXTitle("M[GeV/c^{2}]");
1276 hMassPt->Draw();
1277
595cc7e2 1278 return;
d486095a 1279}
7d9ca2b5 1280//_________________________________________________________________________________________________
1281Int_t AliAnalysisTaskSEDplus::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
1282 //
1283 // checking whether the mother of the particles come from a charm or a bottom quark
1284 //
1285
1286 Int_t pdgGranma = 0;
1287 Int_t mother = 0;
1288 mother = mcPartCandidate->GetMother();
1289 Int_t istep = 0;
1290 Int_t abspdgGranma =0;
1291 Bool_t isFromB=kFALSE;
1292 Bool_t isQuarkFound=kFALSE;
1293 while (mother >0 ){
1294 istep++;
1295 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1296 if (mcGranma){
1297 pdgGranma = mcGranma->GetPdgCode();
1298 abspdgGranma = TMath::Abs(pdgGranma);
1299 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1300 isFromB=kTRUE;
1301 }
1302 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1303 mother = mcGranma->GetMother();
1304 }else{
1305 AliError("Failed casting the mother particle!");
1306 break;
1307 }
1308 }
1309
1310 if(isFromB) return 5;
1311 else return 4;
1312}
af636507 1313//_________________________________________________________________________________________________
1314Float_t AliAnalysisTaskSEDplus::GetTrueImpactParameter(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partDp) const {
1315 // true impact parameter calculation
1316
1317 Double_t vtxTrue[3];
1318 mcHeader->GetVertex(vtxTrue);
1319 Double_t origD[3];
94da94c4 1320 partDp->XvYvZv(origD);
af636507 1321 Short_t charge=partDp->Charge();
1322 Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
1323 for(Int_t iDau=0; iDau<3; iDau++){
1324 pXdauTrue[iDau]=0.;
1325 pYdauTrue[iDau]=0.;
1326 pZdauTrue[iDau]=0.;
1327 }
1328
1329 Int_t nDau=partDp->GetNDaughters();
1330 Int_t labelFirstDau = partDp->GetDaughter(0);
1331 if(nDau==3){
1332 for(Int_t iDau=0; iDau<3; iDau++){
1333 Int_t ind = labelFirstDau+iDau;
1334 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
3de19caa 1335 if(!part){
1336 AliError("Daughter particle not found in MC array");
1337 return 99999.;
1338 }
af636507 1339 pXdauTrue[iDau]=part->Px();
1340 pYdauTrue[iDau]=part->Py();
1341 pZdauTrue[iDau]=part->Pz();
1342 }
1343 }else if(nDau==2){
1344 Int_t theDau=0;
1345 for(Int_t iDau=0; iDau<2; iDau++){
1346 Int_t ind = labelFirstDau+iDau;
1347 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
3de19caa 1348 if(!part){
1349 AliError("Daughter particle not found in MC array");
1350 return 99999.;
1351 }
af636507 1352 Int_t pdgCode=TMath::Abs(part->GetPdgCode());
1353 if(pdgCode==211 || pdgCode==321){
1354 pXdauTrue[theDau]=part->Px();
1355 pYdauTrue[theDau]=part->Py();
1356 pZdauTrue[theDau]=part->Pz();
1357 ++theDau;
1358 }else{
1359 Int_t nDauRes=part->GetNDaughters();
1360 if(nDauRes==2){
1361 Int_t labelFirstDauRes = part->GetDaughter(0);
1362 for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
1363 Int_t indDR = labelFirstDauRes+iDauRes;
1364 AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
3de19caa 1365 if(!partDR){
1366 AliError("Daughter particle not found in MC array");
1367 return 99999.;
1368 }
1369
1370 Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
1371 if(pdgCodeDR==211 || pdgCodeDR==321){
af636507 1372 pXdauTrue[theDau]=partDR->Px();
1373 pYdauTrue[theDau]=partDR->Py();
1374 pZdauTrue[theDau]=partDR->Pz();
1375 ++theDau;
1376 }
1377 }
1378 }
1379 }
1380 }
3de19caa 1381 if(theDau!=3){
1382 AliError("Wrong number of decay prongs");
1383 return 99999.;
1384 }
af636507 1385 }
1386
1387 Double_t d0dummy[3]={0.,0.,0.};
1388 AliAODRecoDecayHF aodDplusMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
1389 return aodDplusMC.ImpParXY();
1390
1391}