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