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