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