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