]> 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
Coverity fixes (Ivana)
[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>
d486095a 33#include <TH1F.h>
4afc48a2 34#include <TDatabasePDG.h>
b557eb43 35
36#include "AliAnalysisManager.h"
4c230f6d 37#include "AliRDHFCutsDplustoKpipi.h"
b557eb43 38#include "AliAODHandler.h"
d486095a 39#include "AliAODEvent.h"
40#include "AliAODVertex.h"
41#include "AliAODTrack.h"
42#include "AliAODMCHeader.h"
43#include "AliAODMCParticle.h"
44#include "AliAODRecoDecayHF3Prong.h"
45#include "AliAnalysisVertexingHF.h"
46#include "AliAnalysisTaskSE.h"
47#include "AliAnalysisTaskSEDplus.h"
a96083b9 48#include "AliNormalizationCounter.h"
d486095a 49
50ClassImp(AliAnalysisTaskSEDplus)
51
52
53//________________________________________________________________________
54AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
55AliAnalysisTaskSE(),
595cc7e2 56 fOutput(0),
57 fHistNEvents(0),
58 fPtVsMass(0),
59 fPtVsMassTC(0),
60 fYVsPt(0),
61 fYVsPtTC(0),
62 fYVsPtSig(0),
63 fYVsPtSigTC(0),
64 fNtupleDplus(0),
65 fUpmasslimit(1.965),
66 fLowmasslimit(1.765),
67 fNPtBins(0),
68 fBinWidth(0.002),
69 fListCuts(0),
70 fRDCutsProduction(0),
71 fRDCutsAnalysis(0),
a96083b9 72 fCounter(0),
595cc7e2 73 fFillNtuple(kFALSE),
74 fReadMC(kFALSE),
a96083b9 75 fUseStrangeness(kFALSE),
595cc7e2 76 fDoLS(kFALSE)
d486095a 77{
82bb8d4b 78 // Default constructor
d486095a 79}
80
81//________________________________________________________________________
4c230f6d 82AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,AliRDHFCutsDplustoKpipi *dpluscutsana,AliRDHFCutsDplustoKpipi *dpluscutsprod,Bool_t fillNtuple):
d486095a 83AliAnalysisTaskSE(name),
4afc48a2 84fOutput(0),
85fHistNEvents(0),
595cc7e2 86fPtVsMass(0),
87fPtVsMassTC(0),
88fYVsPt(0),
89fYVsPtTC(0),
90fYVsPtSig(0),
91fYVsPtSigTC(0),
d486095a 92fNtupleDplus(0),
82bb8d4b 93fUpmasslimit(1.965),
94fLowmasslimit(1.765),
95fNPtBins(0),
993bfba1 96fBinWidth(0.002),
4c230f6d 97fListCuts(0),
98fRDCutsProduction(dpluscutsprod),
99fRDCutsAnalysis(dpluscutsana),
a96083b9 100fCounter(0),
1f4e9722 101fFillNtuple(fillNtuple),
4afc48a2 102fReadMC(kFALSE),
a96083b9 103fUseStrangeness(kFALSE),
4c230f6d 104fDoLS(kFALSE)
d486095a 105{
4c230f6d 106 //
107 // Standrd constructor
108 //
109 //Double_t ptlim[5]={0.,2.,3.,5,9999999.};
110 //SetPtBinLimit(5, ptlim);
111 SetPtBinLimit(fRDCutsAnalysis->GetNPtBins()+1,fRDCutsAnalysis->GetPtBinLimits());
d486095a 112 // Default constructor
82bb8d4b 113 // Output slot #1 writes into a TList container
d486095a 114 DefineOutput(1,TList::Class()); //My private output
4c230f6d 115 // Output slot #2 writes cut to private output
116 // DefineOutput(2,AliRDHFCutsDplustoKpipi::Class());
117 DefineOutput(2,TList::Class());
a96083b9 118// Output slot #3 writes cut to private output
119 DefineOutput(3,AliNormalizationCounter::Class());
120
1f4e9722 121 if(fFillNtuple){
a96083b9 122 // Output slot #4 writes into a TNtuple container
123 DefineOutput(4,TNtuple::Class()); //My private output
1f4e9722 124 }
d486095a 125}
126
127//________________________________________________________________________
128AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
129{
4c230f6d 130 //
d486095a 131 // Destructor
4c230f6d 132 //
d486095a 133 if (fOutput) {
134 delete fOutput;
135 fOutput = 0;
136 }
137c7bcd 137 if(fHistNEvents){
138 delete fHistNEvents;
139 fHistNEvents=0;
140 }
37fc80e6 141 for(Int_t i=0;i<3;i++){
142 if(fHistCentrality[i]){delete fHistCentrality[i]; fHistCentrality[i]=0;}
143 }
137c7bcd 144 for(Int_t i=0;i<3*fNPtBins;i++){
145 if(fMassHist[i]){ delete fMassHist[i]; fMassHist[i]=0;}
146 if(fCosPHist[i]){ delete fCosPHist[i]; fCosPHist[i]=0;}
147 if(fDLenHist[i]){ delete fDLenHist[i]; fDLenHist[i]=0;}
148 if(fSumd02Hist[i]){ delete fSumd02Hist[i]; fSumd02Hist[i]=0;}
149 if(fSigVertHist[i]){ delete fSigVertHist[i]; fSigVertHist[i]=0;}
150 if(fPtMaxHist[i]){ delete fPtMaxHist[i]; fPtMaxHist[i]=0;}
151 if(fDCAHist[i]){ delete fDCAHist[i]; fDCAHist[i]=0;}
152 if(fMassHistTC[i]){ delete fMassHistTC[i]; fMassHistTC[i]=0;}
bb69f127 153 if(fMassHistTCPlus[i]){ delete fMassHistTCPlus[i]; fMassHistTCPlus[i]=0;}
154 if(fMassHistTCMinus[i]){ delete fMassHistTCMinus[i]; fMassHistTCMinus[i]=0;}
4c230f6d 155
137c7bcd 156 if(fMassHistLS[i]){ delete fMassHistLS[i]; fMassHistLS[i]=0;}
157 if(fCosPHistLS[i]){ delete fCosPHistLS[i]; fCosPHistLS[i]=0;}
158 if(fDLenHistLS[i]){ delete fDLenHistLS[i]; fDLenHistLS[i]=0;}
159 if(fSumd02HistLS[i]){ delete fSumd02HistLS[i]; fSumd02HistLS[i]=0;}
160 if(fSigVertHistLS[i]){ delete fSigVertHistLS[i]; fSigVertHistLS[i]=0;}
161 if(fPtMaxHistLS[i]){ delete fPtMaxHistLS[i]; fPtMaxHistLS[i]=0;}
162 if(fDCAHistLS[i]){ delete fDCAHistLS[i]; fDCAHistLS[i]=0;}
163 if(fMassHistLSTC[i]){ delete fMassHistLSTC[i]; fMassHistLSTC[i]=0;}
164 }
165 if(fPtVsMass){
166 delete fPtVsMass;
167 fPtVsMass=0;
168 }
169 if(fPtVsMassTC){
170 delete fPtVsMassTC;
171 fPtVsMassTC=0;
172 }
173 if(fYVsPt){
174 delete fYVsPt;
175 fYVsPt=0;
176 }
177 if(fYVsPtTC){
178 delete fYVsPtTC;
179 fYVsPtTC=0;
180 }
181 if(fYVsPtSig){
182 delete fYVsPtSig;
183 fYVsPtSig=0;
184 }
185 if(fYVsPtSigTC){
186 delete fYVsPtSigTC;
187 fYVsPtSigTC=0;
188 }
189 if(fNtupleDplus){
190 delete fNtupleDplus;
191 fNtupleDplus=0;
192 }
4c230f6d 193 if (fListCuts) {
194 delete fListCuts;
195 fListCuts = 0;
d486095a 196 }
4c230f6d 197 if(fRDCutsProduction){
198 delete fRDCutsProduction;
199 fRDCutsProduction = 0;
200 }
137c7bcd 201 if(fRDCutsAnalysis){
4c230f6d 202 delete fRDCutsAnalysis;
203 fRDCutsAnalysis = 0;
204 }
205
a96083b9 206 if(fCounter){
207 delete fCounter;
208 fCounter = 0;
209 }
137c7bcd 210
211
d486095a 212}
82bb8d4b 213//_________________________________________________________________
214void AliAnalysisTaskSEDplus::SetMassLimits(Float_t range){
4c230f6d 215 // set invariant mass limits
993bfba1 216 Float_t bw=GetBinWidth();
82bb8d4b 217 fUpmasslimit = 1.865+range;
218 fLowmasslimit = 1.865-range;
993bfba1 219 SetBinWidth(bw);
82bb8d4b 220}
221//_________________________________________________________________
222void AliAnalysisTaskSEDplus::SetMassLimits(Float_t lowlimit, Float_t uplimit){
4c230f6d 223 // set invariant mass limits
82bb8d4b 224 if(uplimit>lowlimit)
225 {
993bfba1 226 Float_t bw=GetBinWidth();
82bb8d4b 227 fUpmasslimit = lowlimit;
228 fLowmasslimit = uplimit;
993bfba1 229 SetBinWidth(bw);
82bb8d4b 230 }
231}
82bb8d4b 232//________________________________________________________________________
4c230f6d 233void AliAnalysisTaskSEDplus::SetPtBinLimit(Int_t n, Float_t* lim){
82bb8d4b 234 // define pt bins for analysis
235 if(n>kMaxPtBins){
236 printf("Max. number of Pt bins = %d\n",kMaxPtBins);
237 fNPtBins=kMaxPtBins;
238 fArrayBinLimits[0]=0.;
239 fArrayBinLimits[1]=2.;
240 fArrayBinLimits[2]=3.;
241 fArrayBinLimits[3]=5.;
242 for(Int_t i=4; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
243 }else{
244 fNPtBins=n-1;
245 fArrayBinLimits[0]=lim[0];
246 for(Int_t i=1; i<fNPtBins+1; i++)
247 if(lim[i]>fArrayBinLimits[i-1]){
248 fArrayBinLimits[i]=lim[i];
249 }
250 else {
251 fArrayBinLimits[i]=fArrayBinLimits[i-1];
252 }
253 for(Int_t i=fNPtBins; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
254 }
255 if(fDebug > 1){
256 printf("Number of Pt bins = %d\n",fNPtBins);
e11ae259 257 for(Int_t i=0; i<fNPtBins; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fArrayBinLimits[i],fArrayBinLimits[i+1]);
82bb8d4b 258 }
259}
993bfba1 260//________________________________________________________________
261void AliAnalysisTaskSEDplus::SetBinWidth(Float_t w){
262 Float_t width=w;
263 Int_t nbins=(Int_t)((fUpmasslimit-fLowmasslimit)/width+0.5);
264 Int_t missingbins=4-nbins%4;
265 nbins=nbins+missingbins;
266 width=(fUpmasslimit-fLowmasslimit)/nbins;
267 if(missingbins!=0){
268 printf("AliAnalysisTaskSEDplus::SetBinWidth: W-bin width of %f will produce histograms not rebinnable by 4. New width set to %f\n",w,width);
269 }
270 else{
271 if(fDebug>1) printf("AliAnalysisTaskSEDplus::SetBinWidth: width set to %f\n",width);
272 }
273 fBinWidth=width;
274}
82bb8d4b 275//_________________________________________________________________
276Double_t AliAnalysisTaskSEDplus::GetPtBinLimit(Int_t ibin){
4c230f6d 277 // get pt bin limit
82bb8d4b 278 if(ibin>fNPtBins)return -1;
279 return fArrayBinLimits[ibin];
280}
993bfba1 281//_________________________________________________________________
282Int_t AliAnalysisTaskSEDplus::GetNBinsHistos(){
283 return (Int_t)((fUpmasslimit-fLowmasslimit)/fBinWidth+0.5);
284}
82bb8d4b 285//_________________________________________________________________
286void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS){
4c230f6d 287 //
288 //
289 // Fill the Like Sign histograms
290 //
82bb8d4b 291
292 //count pos/neg tracks
293 Int_t nPosTrks=0,nNegTrks=0;
294 //counter for particles passing single particle cuts
295 Int_t nspcplus=0;
296 Int_t nspcminus=0;
297
298 for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
299 AliAODTrack *track = aod->GetTrack(it);
300 if(track->Charge()>0){
301 nPosTrks++;
302 if(track->Pt()>=0.4){
303 nspcplus++;
304 }
305 }
306 if(track->Charge()<0)
307 {
308 nNegTrks++;
309 if(track->Pt()>=0.4){
310 nspcminus++;
311 }
312 }
313 }
314
315 Int_t nOStriplets = arrayOppositeSign->GetEntriesFast();
316
317 Int_t nDplusLS=0;
318 Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
319 Int_t index;
320
321 for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
322 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
323 Bool_t unsetvtx=kFALSE;
324 if(!d->GetOwnPrimaryVtx()) {
325 d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
326 unsetvtx=kTRUE;
327 }
1275bf81 328 if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod))nDplusLS++;
82bb8d4b 329 if(unsetvtx) d->UnsetOwnPrimaryVtx();
330 }
331
332 Float_t wei2=0;
333 if(nLikeSign!=0)wei2 = (Float_t)nOStriplets/(Float_t)nLikeSign;
334 Float_t wei3=0;
335 if(nDplusLS!=0)wei3 = (Float_t)nDplusOS/(Float_t)nDplusLS;
336
337 // loop over like sign candidates
338 for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
339 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
340 Bool_t unsetvtx=kFALSE;
341 if(!d->GetOwnPrimaryVtx()) {
342 d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
343 unsetvtx=kTRUE;
344 }
345
1275bf81 346 if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)){
82bb8d4b 347
348 //set tight cuts values
349 Int_t iPtBin=-1;
350 Double_t ptCand = d->Pt();
82bb8d4b 351 for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
352 if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
353 }
354
355 if(iPtBin<0){
356 return;
357 }
358
1275bf81 359 Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
82bb8d4b 360
361 Int_t sign= d->GetCharge();
362 Float_t wei=1;
363 Float_t wei4=1;
364 if(sign>0&&nPosTrks>2&&nspcplus>2) { //wei* should be automatically protected, since to get a triplet there must be at least 3 good tracks in the event
365
366 wei=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.);
367 wei4=3.*(Float_t)nspcminus/((Float_t)nspcplus-2.);
368 }
369
370 if(sign<0&&nNegTrks>2&&nspcminus>2){
371 wei=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.);
372 wei4=3.*(Float_t)nspcplus/((Float_t)nspcminus-2.);
373
374 }
375
376 Float_t invMass = d->InvMassDplus();
ba9ae5b2 377 Double_t dlen=d->DecayLength();
378 Double_t cosp=d->CosPointingAngle();
379 Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
380 Double_t dca=d->GetDCA();
0db12536 381 Double_t sigvert=d->GetSigmaVert();
ba9ae5b2 382 Double_t ptmax=0;
383 for(Int_t i=0;i<3;i++){
384 if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
385 }
386
82bb8d4b 387 index=GetLSHistoIndex(iPtBin);
388 fMassHistLS[index]->Fill(invMass,wei);
389 fMassHistLS[index+1]->Fill(invMass);
390 fMassHistLS[index+2]->Fill(invMass,wei2);
391 fMassHistLS[index+3]->Fill(invMass,wei3);
392 fMassHistLS[index+4]->Fill(invMass,wei4);
ba9ae5b2 393
394 Int_t indexcut=GetHistoIndex(iPtBin);
395 fCosPHistLS[indexcut]->Fill(cosp);
396 fDLenHistLS[indexcut]->Fill(dlen);
397 fSumd02HistLS[indexcut]->Fill(sumD02);
0db12536 398 fSigVertHistLS[indexcut]->Fill(sigvert);
ba9ae5b2 399 fPtMaxHistLS[indexcut]->Fill(ptmax);
400 fDCAHistLS[indexcut]->Fill(dca);
82bb8d4b 401
595cc7e2 402 if(passTightCuts){
82bb8d4b 403 fMassHistLSTC[index]->Fill(invMass,wei);
404 fMassHistLSTC[index+1]->Fill(invMass);
405 fMassHistLSTC[index+2]->Fill(invMass,wei2);
406 fMassHistLSTC[index+3]->Fill(invMass,wei3);
407 fMassHistLSTC[index+4]->Fill(invMass,wei4);
408 }
409 }
410 if(unsetvtx) d->UnsetOwnPrimaryVtx();
411 }
412
413 //printf("------------ N. of positive tracks in Event ----- %d \n", nPosTrks);
414 //printf("------------ N. of negative tracks in Event ----- %d \n", nNegTrks);
415
416 // printf("LS analysis...done\n");
417
418}
d486095a 419
d486095a 420
4c230f6d 421//__________________________________________
422void AliAnalysisTaskSEDplus::Init(){
423 //
424 // Initialization
425 //
d486095a 426 if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
4c230f6d 427
428 //PostData(2,fRDCutsloose);//we should then put those cuts in a tlist if we have more than 1
429 fListCuts=new TList();
90993728 430 AliRDHFCutsDplustoKpipi *production = new AliRDHFCutsDplustoKpipi(*fRDCutsProduction);
431 production->SetName("ProductionCuts");
432 AliRDHFCutsDplustoKpipi *analysis = new AliRDHFCutsDplustoKpipi(*fRDCutsAnalysis);
433 analysis->SetName("AnalysisCuts");
4c230f6d 434
435 fListCuts->Add(production);
436 fListCuts->Add(analysis);
437 PostData(2,fListCuts);
438
d486095a 439 return;
440}
441
442//________________________________________________________________________
443void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
444{
445 // Create the output container
446 //
447 if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");
448
449 // Several histograms are more conveniently managed in a TList
450 fOutput = new TList();
451 fOutput->SetOwner();
1f4e9722 452 fOutput->SetName("OutputHistos");
453
1f4e9722 454 TString hisname;
82bb8d4b 455 Int_t index=0;
456 Int_t indexLS=0;
993bfba1 457 Int_t nbins=GetNBinsHistos();
37fc80e6 458 fHistCentrality[0]=new TH1F("centrality","centrality",100,0.5,30000.5);
459 fHistCentrality[1]=new TH1F("centrality(selectedCent)","centrality(selectedCent)",100,0.5,30000.5);
460 fHistCentrality[2]=new TH1F("centrality(OutofCent)","centrality(OutofCent)",100,0.5,30000.5);
461 for(Int_t i=0;i<3;i++){
462 fHistCentrality[i]->Sumw2();
463 fOutput->Add(fHistCentrality[i]);
464 }
82bb8d4b 465 for(Int_t i=0;i<fNPtBins;i++){
1f4e9722 466
82bb8d4b 467 index=GetHistoIndex(i);
468 indexLS=GetLSHistoIndex(i);
ba9ae5b2 469
82bb8d4b 470 hisname.Form("hMassPt%d",i);
993bfba1 471 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 472 fMassHist[index]->Sumw2();
ba9ae5b2 473 hisname.Form("hCosPAllPt%d",i);
993bfba1 474 fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
ba9ae5b2 475 fCosPHist[index]->Sumw2();
476 hisname.Form("hDLenAllPt%d",i);
993bfba1 477 fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
ba9ae5b2 478 fDLenHist[index]->Sumw2();
479 hisname.Form("hSumd02AllPt%d",i);
993bfba1 480 fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
ba9ae5b2 481 fSumd02Hist[index]->Sumw2();
482 hisname.Form("hSigVertAllPt%d",i);
993bfba1 483 fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
ba9ae5b2 484 fSigVertHist[index]->Sumw2();
485 hisname.Form("hPtMaxAllPt%d",i);
993bfba1 486 fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
ba9ae5b2 487 fPtMaxHist[index]->Sumw2();
488
489 hisname.Form("hDCAAllPt%d",i);
993bfba1 490 fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
ba9ae5b2 491 fDCAHist[index]->Sumw2();
492
493
494
1f4e9722 495 hisname.Form("hMassPt%dTC",i);
993bfba1 496 fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 497 fMassHistTC[index]->Sumw2();
bb69f127 498 hisname.Form("hMassPt%dTCPlus",i);
499 fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
500 fMassHistTCPlus[index]->Sumw2();
501 hisname.Form("hMassPt%dTCMinus",i);
502 fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
503 fMassHistTCMinus[index]->Sumw2();
ba9ae5b2 504
505
506
507
508 hisname.Form("hCosPAllPt%dLS",i);
993bfba1 509 fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
ba9ae5b2 510 fCosPHistLS[index]->Sumw2();
511 hisname.Form("hDLenAllPt%dLS",i);
993bfba1 512 fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
ba9ae5b2 513 fDLenHistLS[index]->Sumw2();
514 hisname.Form("hSumd02AllPt%dLS",i);
993bfba1 515 fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
ba9ae5b2 516 fSumd02HistLS[index]->Sumw2();
517 hisname.Form("hSigVertAllPt%dLS",i);
993bfba1 518 fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
ba9ae5b2 519 fSigVertHistLS[index]->Sumw2();
520 hisname.Form("hPtMaxAllPt%dLS",i);
993bfba1 521 fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
ba9ae5b2 522 fPtMaxHistLS[index]->Sumw2();
523
524 hisname.Form("hDCAAllPt%dLS",i);
993bfba1 525 fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
ba9ae5b2 526 fDCAHistLS[index]->Sumw2();
527
82bb8d4b 528 hisname.Form("hLSPt%dLC",i);
993bfba1 529 fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 530 fMassHistLS[indexLS]->Sumw2();
531
ba9ae5b2 532 hisname.Form("hLSPt%dTC",i);
993bfba1 533 fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
ba9ae5b2 534 fMassHistLSTC[indexLS]->Sumw2();
535
536
537
82bb8d4b 538 index=GetSignalHistoIndex(i);
539 indexLS++;
540 hisname.Form("hSigPt%d",i);
993bfba1 541 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 542 fMassHist[index]->Sumw2();
ba9ae5b2 543 hisname.Form("hCosPSigPt%d",i);
993bfba1 544 fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
ba9ae5b2 545 fCosPHist[index]->Sumw2();
546 hisname.Form("hDLenSigPt%d",i);
993bfba1 547 fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
ba9ae5b2 548 fDLenHist[index]->Sumw2();
549 hisname.Form("hSumd02SigPt%d",i);
993bfba1 550 fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
ba9ae5b2 551 fSumd02Hist[index]->Sumw2();
552 hisname.Form("hSigVertSigPt%d",i);
993bfba1 553 fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
ba9ae5b2 554 fSigVertHist[index]->Sumw2();
555 hisname.Form("hPtMaxSigPt%d",i);
993bfba1 556 fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
ba9ae5b2 557 fPtMaxHist[index]->Sumw2();
558
559 hisname.Form("hDCASigPt%d",i);
993bfba1 560 fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
ba9ae5b2 561 fDCAHist[index]->Sumw2();
562
563
1f4e9722 564 hisname.Form("hSigPt%dTC",i);
993bfba1 565 fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 566 fMassHistTC[index]->Sumw2();
bb69f127 567 hisname.Form("hSigPt%dTCPlus",i);
568 fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
569 fMassHistTCPlus[index]->Sumw2();
570 hisname.Form("hSigPt%dTCMinus",i);
571 fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
572 fMassHistTCMinus[index]->Sumw2();
ba9ae5b2 573
82bb8d4b 574 hisname.Form("hLSPt%dLCnw",i);
993bfba1 575 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 576 fMassHistLS[indexLS]->Sumw2();
577 hisname.Form("hLSPt%dTCnw",i);
993bfba1 578 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 579 fMassHistLSTC[indexLS]->Sumw2();
580
ba9ae5b2 581
582
583 hisname.Form("hCosPSigPt%dLS",i);
993bfba1 584 fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
ba9ae5b2 585 fCosPHistLS[index]->Sumw2();
586 hisname.Form("hDLenSigPt%dLS",i);
993bfba1 587 fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
ba9ae5b2 588 fDLenHistLS[index]->Sumw2();
589 hisname.Form("hSumd02SigPt%dLS",i);
993bfba1 590 fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
ba9ae5b2 591 fSumd02HistLS[index]->Sumw2();
592 hisname.Form("hSigVertSigPt%dLS",i);
993bfba1 593 fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
ba9ae5b2 594 fSigVertHistLS[index]->Sumw2();
595 hisname.Form("hPtMaxSigPt%dLS",i);
993bfba1 596 fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
ba9ae5b2 597 fPtMaxHistLS[index]->Sumw2();
598
599 hisname.Form("hDCASigPt%dLS",i);
993bfba1 600 fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
ba9ae5b2 601 fDCAHistLS[index]->Sumw2();
602
603
604
82bb8d4b 605 index=GetBackgroundHistoIndex(i);
606 indexLS++;
607 hisname.Form("hBkgPt%d",i);
993bfba1 608 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 609 fMassHist[index]->Sumw2();
ba9ae5b2 610 hisname.Form("hCosPBkgPt%d",i);
993bfba1 611 fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
ba9ae5b2 612 fCosPHist[index]->Sumw2();
613 hisname.Form("hDLenBkgPt%d",i);
993bfba1 614 fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
ba9ae5b2 615 fDLenHist[index]->Sumw2();
616 hisname.Form("hSumd02BkgPt%d",i);
993bfba1 617 fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
ba9ae5b2 618 fSumd02Hist[index]->Sumw2();
619 hisname.Form("hSigVertBkgPt%d",i);
993bfba1 620 fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
ba9ae5b2 621 fSigVertHist[index]->Sumw2();
622 hisname.Form("hPtMaxBkgPt%d",i);
993bfba1 623 fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
ba9ae5b2 624 fPtMaxHist[index]->Sumw2();
625
626 hisname.Form("hDCABkgPt%d",i);
993bfba1 627 fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
ba9ae5b2 628 fDCAHist[index]->Sumw2();
629
630
1f4e9722 631 hisname.Form("hBkgPt%dTC",i);
993bfba1 632 fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 633 fMassHistTC[index]->Sumw2();
bb69f127 634 hisname.Form("hBkgPt%dTCPlus",i);
635 fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
636 fMassHistTCPlus[index]->Sumw2();
637 hisname.Form("hBkgPt%dTCMinus",i);
638 fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
639 fMassHistTCMinus[index]->Sumw2();
ba9ae5b2 640
82bb8d4b 641 hisname.Form("hLSPt%dLCntrip",i);
993bfba1 642 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 643 fMassHistLS[indexLS]->Sumw2();
644 hisname.Form("hLSPt%dTCntrip",i);
993bfba1 645 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 646 fMassHistLSTC[indexLS]->Sumw2();
647
ba9ae5b2 648
649 hisname.Form("hCosPBkgPt%dLS",i);
993bfba1 650 fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
ba9ae5b2 651 fCosPHistLS[index]->Sumw2();
652 hisname.Form("hDLenBkgPt%dLS",i);
993bfba1 653 fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
ba9ae5b2 654 fDLenHistLS[index]->Sumw2();
655 hisname.Form("hSumd02BkgPt%dLS",i);
993bfba1 656 fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
ba9ae5b2 657 fSumd02HistLS[index]->Sumw2();
658 hisname.Form("hSigVertBkgPt%dLS",i);
993bfba1 659 fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
ba9ae5b2 660 fSigVertHistLS[index]->Sumw2();
661 hisname.Form("hPtMaxBkgPt%dLS",i);
993bfba1 662 fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
ba9ae5b2 663 fPtMaxHistLS[index]->Sumw2();
664 hisname.Form("hDCABkgPt%dLS",i);
993bfba1 665 fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
ba9ae5b2 666 fDCAHistLS[index]->Sumw2();
667
668
82bb8d4b 669 indexLS++;
670 hisname.Form("hLSPt%dLCntripsinglecut",i);
993bfba1 671 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 672 fMassHistLS[indexLS]->Sumw2();
673 hisname.Form("hLSPt%dTCntripsinglecut",i);
993bfba1 674 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 675 fMassHistLSTC[indexLS]->Sumw2();
676
677 indexLS++;
678 hisname.Form("hLSPt%dLCspc",i);
993bfba1 679 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 680 fMassHistLS[indexLS]->Sumw2();
681 hisname.Form("hLSPt%dTCspc",i);
993bfba1 682 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
82bb8d4b 683 fMassHistLSTC[indexLS]->Sumw2();
1f4e9722 684 }
8c34bd86 685
82bb8d4b 686 for(Int_t i=0; i<3*fNPtBins; i++){
687 fOutput->Add(fMassHist[i]);
ba9ae5b2 688 fOutput->Add(fCosPHist[i]);
689 fOutput->Add(fDLenHist[i]);
690 fOutput->Add(fSumd02Hist[i]);
691 fOutput->Add(fSigVertHist[i]);
692 fOutput->Add(fPtMaxHist[i]);
693 fOutput->Add(fDCAHist[i]);
82bb8d4b 694 fOutput->Add(fMassHistTC[i]);
bb69f127 695 fOutput->Add(fMassHistTCPlus[i]);
696 fOutput->Add(fMassHistTCMinus[i]);
82bb8d4b 697 }
ba9ae5b2 698 for(Int_t i=0; i<3*fNPtBins&&fDoLS; i++){
699 fOutput->Add(fCosPHistLS[i]);
700 fOutput->Add(fDLenHistLS[i]);
701 fOutput->Add(fSumd02HistLS[i]);
702 fOutput->Add(fSigVertHistLS[i]);
703 fOutput->Add(fPtMaxHistLS[i]);
704 fOutput->Add(fDCAHistLS[i]);
705 }
82bb8d4b 706 for(Int_t i=0; i<5*fNPtBins&&fDoLS; i++){
707 fOutput->Add(fMassHistLS[i]);
708 fOutput->Add(fMassHistLSTC[i]);
709 }
ba9ae5b2 710
26da4c42 711
37fc80e6 712 fHistNEvents = new TH1F("fHistNEvents", "number of events ",8,-0.5,7.5);
de1b82ff 713 fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
26da4c42 714 fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with good vertex");
de1b82ff 715 fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents with PbPb HM trigger");
716 fHistNEvents->GetXaxis()->SetBinLabel(4,"no. of Rejected pileup events");
717 fHistNEvents->GetXaxis()->SetBinLabel(5,"no. of candidate");
718 fHistNEvents->GetXaxis()->SetBinLabel(6,"no. of D+ after loose cuts");
719 fHistNEvents->GetXaxis()->SetBinLabel(7,"no. of D+ after tight cuts");
37fc80e6 720 fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events");
721
26da4c42 722 fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
723
ba9ae5b2 724 fHistNEvents->Sumw2();
725 fHistNEvents->SetMinimum(0);
726 fOutput->Add(fHistNEvents);
595cc7e2 727
4d797793 728 fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
729 fPtVsMassTC=new TH2F("hPtVsMassTC","PtVsMass (analysis cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
595cc7e2 730 fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
731 fYVsPtTC=new TH2F("hYVsPtTC","YvsPt (analysis cuts)",40,0.,20.,80,-2.,2.);
732 fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
733 fYVsPtSigTC=new TH2F("hYVsPtSigTC","YvsPt (MC, only Sig, analysis cuts)",40,0.,20.,80,-2.,2.);
734
735 fOutput->Add(fPtVsMass);
736 fOutput->Add(fPtVsMassTC);
737 fOutput->Add(fYVsPt);
738 fOutput->Add(fYVsPtTC);
739 fOutput->Add(fYVsPtSig);
740 fOutput->Add(fYVsPtSigTC);
4afc48a2 741
a96083b9 742
743 //Counter for Normalization
e0fdfa52 744 TString normName="NormalizationCounter";
745 AliAnalysisDataContainer *cont = GetOutputSlot(3)->GetContainer();
746 if(cont)normName=(TString)cont->GetName();
747 fCounter = new AliNormalizationCounter(normName.Data());
de1b82ff 748 fCounter->SetRejectPileUp(fRDCutsProduction->GetOptPileUp());
a96083b9 749
1f4e9722 750 if(fFillNtuple){
e0fdfa52 751 OpenFile(4); // 4 is the slot number of the ntuple
ba9ae5b2 752
753 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");
754
1f4e9722 755 }
ba9ae5b2 756
d486095a 757 return;
758}
759
760//________________________________________________________________________
761void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
762{
763 // Execute analysis for current event:
764 // heavy flavor candidates association to MC truth
765
4c230f6d 766 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
4afc48a2 767
a96083b9 768
769
b557eb43 770 TClonesArray *array3Prong = 0;
4afc48a2 771 TClonesArray *arrayLikeSign =0;
b557eb43 772 if(!aod && AODEvent() && IsStandardAOD()) {
773 // In case there is an AOD handler writing a standard AOD, use the AOD
774 // event in memory rather than the input (ESD) event.
775 aod = dynamic_cast<AliAODEvent*> (AODEvent());
4c230f6d 776 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
777 // have to taken from the AOD event hold by the AliAODExtension
b557eb43 778 AliAODHandler* aodHandler = (AliAODHandler*)
779 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
780 if(aodHandler->GetExtensions()) {
781 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
782 AliAODEvent *aodFromExt = ext->GetAOD();
783 array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
4afc48a2 784 arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
b557eb43 785 }
dc222f77 786 } else if(aod) {
b557eb43 787 array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
4afc48a2 788 arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
b557eb43 789 }
8931c313 790
90993728 791 if(!aod || !array3Prong) {
d486095a 792 printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
793 return;
794 }
3a3c90f4 795 if(!arrayLikeSign && fDoLS) {
4afc48a2 796 printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
82bb8d4b 797 return;
4afc48a2 798 }
799
7c23877d 800
801 // fix for temporary bug in ESDfilter
802 // the AODs with null vertex pointer didn't pass the PhysSel
a96083b9 803 if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
804 fCounter->StoreEvent(aod,fReadMC);
de1b82ff 805 fHistNEvents->Fill(0); // count event
37fc80e6 806
807 Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
808 Float_t centrality=aod->GetNTracks();//fRDCutsAnalysis->GetCentrality(aod);
809 fHistCentrality[0]->Fill(centrality);
de1b82ff 810 // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
811 TString trigclass=aod->GetFiredTriggerClasses();
812 if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")||trigclass.Contains("C0SMH-B-NOPF-ALL")) fHistNEvents->Fill(2);
813 if(fRDCutsAnalysis->GetWhyRejection()==1)fHistNEvents->Fill(3);
37fc80e6 814 if(fRDCutsAnalysis->GetWhyRejection()==2){fHistNEvents->Fill(7);fHistCentrality[2]->Fill(centrality);}
815
816 // Post the data already here
de1b82ff 817 PostData(1,fOutput);
37fc80e6 818 if(!isEvSel)return;
819
820 fHistCentrality[1]->Fill(centrality);
821 fHistNEvents->Fill(1);
822
4afc48a2 823 TClonesArray *arrayMC=0;
824 AliAODMCHeader *mcHeader=0;
d486095a 825
826 // AOD primary vertex
1f4e9722 827 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
828 // vtx1->Print();
26da4c42 829 TString primTitle = vtx1->GetTitle();
830 //if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0)fHistNEvents->Fill(2);
831
d486095a 832 // load MC particles
4afc48a2 833 if(fReadMC){
834
835 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
836 if(!arrayMC) {
837 printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
90993728 838 return;
4afc48a2 839 }
1f4e9722 840
d486095a 841 // load MC header
4afc48a2 842 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
843 if(!mcHeader) {
82bb8d4b 844 printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
845 return;
4afc48a2 846 }
d486095a 847 }
4afc48a2 848
1f4e9722 849 Int_t n3Prong = array3Prong->GetEntriesFast();
850 if(fDebug>1) printf("Number of D+->Kpipi: %d\n",n3Prong);
4afc48a2 851
852
853 Int_t nOS=0;
82bb8d4b 854 Int_t index;
1f4e9722 855 Int_t pdgDgDplustoKpipi[3]={321,211,211};
4c230f6d 856 // 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
857 //Double_t *cutsDplus = new (Double_t*)fRDCuts->GetCuts();
a96083b9 858 Int_t nSelectedloose=0,nSelectedtight=0;
1f4e9722 859 for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
860 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
de1b82ff 861 fHistNEvents->Fill(4);
1f4e9722 862
863 Bool_t unsetvtx=kFALSE;
864 if(!d->GetOwnPrimaryVtx()){
865 d->SetOwnPrimaryVtx(vtx1);
866 unsetvtx=kTRUE;
867 }
4afc48a2 868
1275bf81 869 if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)) {
a96083b9 870
26da4c42 871
a96083b9 872
82bb8d4b 873 Int_t iPtBin = -1;
1f4e9722 874 Double_t ptCand = d->Pt();
4c230f6d 875
82bb8d4b 876 for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
877 if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
1f4e9722 878 }
82bb8d4b 879
1275bf81 880 Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
a96083b9 881
26da4c42 882
4afc48a2 883 Int_t labDp=-1;
884 Float_t deltaPx=0.;
885 Float_t deltaPy=0.;
886 Float_t deltaPz=0.;
887 Float_t truePt=0.;
888 Float_t xDecay=0.;
889 Float_t yDecay=0.;
890 Float_t zDecay=0.;
891 Float_t pdgCode=-2;
892 if(fReadMC){
893 labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
894 if(labDp>=0){
895 AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
896 AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
897 deltaPx=partDp->Px()-d->Px();
898 deltaPy=partDp->Py()-d->Py();
899 deltaPz=partDp->Pz()-d->Pz();
900 truePt=partDp->Pt();
901 xDecay=dg0->Xv();
902 yDecay=dg0->Yv();
903 zDecay=dg0->Zv();
904 pdgCode=TMath::Abs(partDp->GetPdgCode());
905 }else{
906 pdgCode=-1;
907 }
908 }
1f4e9722 909 Double_t invMass=d->InvMassDplus();
595cc7e2 910 Double_t rapid=d->YDplus();
911 fYVsPt->Fill(ptCand,rapid);
912 if(passTightCuts) fYVsPtTC->Fill(ptCand,rapid);
913 Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
914 if(isFidAcc){
915 fPtVsMass->Fill(invMass,ptCand);
916 if(passTightCuts) fPtVsMassTC->Fill(invMass,ptCand);
917 }
ba9ae5b2 918 Float_t tmp[24];
4afc48a2 919 if(fFillNtuple){
920 tmp[0]=pdgCode;
921 tmp[1]=deltaPx;
922 tmp[2]=deltaPy;
923 tmp[3]=deltaPz;
924 tmp[4]=truePt;
925 tmp[5]=xDecay;
926 tmp[6]=yDecay;
927 tmp[7]=zDecay;
928 tmp[8]=d->PtProng(0);
929 tmp[9]=d->PtProng(1);
930 tmp[10]=d->PtProng(2);
931 tmp[11]=d->Pt();
932 tmp[12]=d->CosPointingAngle();
933 tmp[13]=d->DecayLength();
934 tmp[14]=d->Xv();
935 tmp[15]=d->Yv();
936 tmp[16]=d->Zv();
937 tmp[17]=d->InvMassDplus();
938 tmp[18]=d->GetSigmaVert();
939 tmp[19]=d->Getd0Prong(0);
940 tmp[20]=d->Getd0Prong(1);
ba9ae5b2 941 tmp[21]=d->Getd0Prong(2);
942 tmp[22]=d->GetDCA();
943 tmp[23]=d->Prodd0d0();
4afc48a2 944 fNtupleDplus->Fill(tmp);
de1b82ff 945 PostData(4,fNtupleDplus);
4afc48a2 946 }
ba9ae5b2 947 Double_t dlen=d->DecayLength();
948 Double_t cosp=d->CosPointingAngle();
949 Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
0db12536 950 Double_t dca=d->GetDCA();
951 Double_t sigvert=d->GetSigmaVert();
952 Double_t ptmax=0;
ba9ae5b2 953 for(Int_t i=0;i<3;i++){
954 if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
955 }
82bb8d4b 956 if(iPtBin>=0){
957
958 index=GetHistoIndex(iPtBin);
595cc7e2 959 if(isFidAcc){
de1b82ff 960 fHistNEvents->Fill(5);
a96083b9 961 nSelectedloose++;
595cc7e2 962 fMassHist[index]->Fill(invMass);
963 fCosPHist[index]->Fill(cosp);
964 fDLenHist[index]->Fill(dlen);
965 fSumd02Hist[index]->Fill(sumD02);
966 fSigVertHist[index]->Fill(sigvert);
967 fPtMaxHist[index]->Fill(ptmax);
968 fDCAHist[index]->Fill(dca);
969
de1b82ff 970 if(passTightCuts){ fHistNEvents->Fill(6);
a96083b9 971 nSelectedtight++;
595cc7e2 972 fMassHistTC[index]->Fill(invMass);
bb69f127 973 if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
974 else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
595cc7e2 975 }
82bb8d4b 976 }
977
978 if(fReadMC){
979 if(labDp>=0) {
980 index=GetSignalHistoIndex(iPtBin);
595cc7e2 981 if(isFidAcc){
a96083b9 982 Float_t factor[3]={1.,1.,1.};
983 if(fUseStrangeness){
984 for(Int_t iprong=0;iprong<3;iprong++){
985 AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
986 Int_t labd= trad->GetLabel();
987 if(labd>=0){
988 AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
989 if(dau){
990 Int_t labm = dau->GetMother();
991 if(labm>=0){
992 AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
993 if(mot){
994 if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
995 if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
996 else factor[iprong]=1./.6;
997 // fNentries->Fill(12);
998 }
999 if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
1000 factor[iprong]=1./0.25;
1001 // fNentries->Fill(13);
1002 }//if 3122
1003 }//if(mot)
1004 }//if labm>0
1005 }//if(dau)
1006 }//if labd>=0
1007 }//prong loop
1008 }
1009 Float_t fact=1.;for(Int_t k=0;k<3;k++)fact=fact*factor[k];
595cc7e2 1010 fMassHist[index]->Fill(invMass);
a96083b9 1011 fCosPHist[index]->Fill(cosp,fact);
1012 fDLenHist[index]->Fill(dlen,fact);
1013 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];
1014 fSumd02Hist[index]->Fill(sumd02s);
1015 fSigVertHist[index]->Fill(sigvert,fact);
1016 fPtMaxHist[index]->Fill(ptmax,fact);
1017 fDCAHist[index]->Fill(dca,fact);
595cc7e2 1018 if(passTightCuts){
1019 fMassHistTC[index]->Fill(invMass);
bb69f127 1020 if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
1021 else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
595cc7e2 1022 }
1023 }
1024 fYVsPtSig->Fill(ptCand,rapid);
1025 if(passTightCuts) fYVsPtSigTC->Fill(ptCand,rapid);
82bb8d4b 1026 }else{
1027 index=GetBackgroundHistoIndex(iPtBin);
595cc7e2 1028 if(isFidAcc){
a96083b9 1029 Float_t factor[3]={1.,1.,1.};
1030 if(fUseStrangeness){
1031 for(Int_t iprong=0;iprong<3;iprong++){
1032 AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
1033 Int_t labd= trad->GetLabel();
1034 if(labd>=0){
1035 AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
1036 if(dau){
1037 Int_t labm = dau->GetMother();
1038 if(labm>=0){
1039 AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
1040 if(mot){
1041 if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
1042 if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
1043 else factor[iprong]=1./.6;
1044 // fNentries->Fill(12);
1045 }
1046 if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
1047 factor[iprong]=1./0.25;
1048 // fNentries->Fill(13);
1049 }//if 3122
1050 }//if(mot)
1051 }//if labm>0
1052 }//if(dau)
1053 }//if labd>=0
1054 }//prong loop
1055 }
1056 Float_t fact=1.;for(Int_t k=0;k<3;k++)fact=fact*factor[k];
595cc7e2 1057 fMassHist[index]->Fill(invMass);
a96083b9 1058 fCosPHist[index]->Fill(cosp,fact);
1059 fDLenHist[index]->Fill(dlen,fact);
1060 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];
1061 fSumd02Hist[index]->Fill(sumd02s);
1062 fSigVertHist[index]->Fill(sigvert,fact);
1063 fPtMaxHist[index]->Fill(ptmax,fact);
1064 fDCAHist[index]->Fill(dca,fact);
595cc7e2 1065 if(passTightCuts){
1066 fMassHistTC[index]->Fill(invMass);
bb69f127 1067 if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
1068 else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
595cc7e2 1069 }
82bb8d4b 1070 }
4afc48a2 1071 }
fc8d975b 1072 }
a96083b9 1073 }
d486095a 1074 }
1f4e9722 1075 if(unsetvtx) d->UnsetOwnPrimaryVtx();
1076 }
a96083b9 1077 fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
1078 fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
4c230f6d 1079
4afc48a2 1080 //start LS analysis
1081 if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
82bb8d4b 1082
4afc48a2 1083 PostData(1,fOutput);
a96083b9 1084 PostData(3,fCounter);
d486095a 1085 return;
1086}
1087
4afc48a2 1088
1089
82bb8d4b 1090//________________________________________________________________________
d486095a 1091void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
1092{
1093 // Terminate analysis
1094 //
1095 if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
82bb8d4b 1096
d486095a 1097 fOutput = dynamic_cast<TList*> (GetOutputData(1));
1098 if (!fOutput) {
1099 printf("ERROR: fOutput not available\n");
1100 return;
1101 }
595cc7e2 1102 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
82bb8d4b 1103
595cc7e2 1104 TString hisname;
1105 Int_t index=0;
ba9ae5b2 1106
595cc7e2 1107 for(Int_t i=0;i<fNPtBins;i++){
82bb8d4b 1108 index=GetHistoIndex(i);
ba9ae5b2 1109
37fc80e6 1110 hisname.Form("hMassPt%dTC",i);
82bb8d4b 1111 fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
37fc80e6 1112 }
ba9ae5b2 1113
13808a30 1114 TCanvas *c1=new TCanvas("c1","D+ invariant mass distribution",500,500);
1115 c1->cd();
1116 TH1F *hMassPt=(TH1F*)fOutput->FindObject("hMassPt3TC");
1117 hMassPt->SetLineColor(kBlue);
1118 hMassPt->SetXTitle("M[GeV/c^{2}]");
1119 hMassPt->Draw();
1120
595cc7e2 1121 return;
d486095a 1122}