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