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