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