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