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