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