]> 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
Fix warnings
[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
16/////////////////////////////////////////////////////////////
17//
18// AliAnalysisTaskSE for the extraction of signal(e.g D+) of heavy flavor
19// decay candidates with the MC truth.
4afc48a2 20// Authors: Renu Bala, bala@to.infn.it
21// F. Prino, prino@to.infn.it
22// G. Ortona, ortona@to.infn.it
d486095a 23/////////////////////////////////////////////////////////////
24
25#include <TClonesArray.h>
26#include <TNtuple.h>
13808a30 27#include <TCanvas.h>
d486095a 28#include <TList.h>
1f4e9722 29#include <TString.h>
d486095a 30#include <TH1F.h>
4afc48a2 31#include <TDatabasePDG.h>
b557eb43 32
33#include "AliAnalysisManager.h"
34#include "AliAODHandler.h"
d486095a 35#include "AliAODEvent.h"
36#include "AliAODVertex.h"
37#include "AliAODTrack.h"
38#include "AliAODMCHeader.h"
39#include "AliAODMCParticle.h"
40#include "AliAODRecoDecayHF3Prong.h"
41#include "AliAnalysisVertexingHF.h"
42#include "AliAnalysisTaskSE.h"
43#include "AliAnalysisTaskSEDplus.h"
44
45ClassImp(AliAnalysisTaskSEDplus)
46
47
48//________________________________________________________________________
49AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
50AliAnalysisTaskSE(),
51fOutput(0),
4afc48a2 52fHistNEvents(0),
d486095a 53fNtupleDplus(0),
82bb8d4b 54fUpmasslimit(1.965),
55fLowmasslimit(1.765),
56fNPtBins(0),
1f4e9722 57fFillNtuple(kFALSE),
4afc48a2 58fReadMC(kFALSE),
59fDoLS(kFALSE),
d486095a 60fVHF(0)
61{
82bb8d4b 62 // Default constructor
d486095a 63}
64
65//________________________________________________________________________
1f4e9722 66AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,Bool_t fillNtuple):
d486095a 67AliAnalysisTaskSE(name),
4afc48a2 68fOutput(0),
69fHistNEvents(0),
d486095a 70fNtupleDplus(0),
82bb8d4b 71fUpmasslimit(1.965),
72fLowmasslimit(1.765),
73fNPtBins(0),
1f4e9722 74fFillNtuple(fillNtuple),
4afc48a2 75fReadMC(kFALSE),
76fDoLS(kFALSE),
d486095a 77fVHF(0)
78{
82bb8d4b 79 Double_t ptlim[5]={0.,2.,3.,5,9999999.};
80 SetPtBinLimit(5, ptlim);
d486095a 81 // Default constructor
82bb8d4b 82 // Output slot #1 writes into a TList container
d486095a 83 DefineOutput(1,TList::Class()); //My private output
1f4e9722 84
85 if(fFillNtuple){
86 // Output slot #2 writes into a TNtuple container
87 DefineOutput(2,TNtuple::Class()); //My private output
1f4e9722 88 }
d486095a 89}
90
91//________________________________________________________________________
92AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
93{
94 // Destructor
95 if (fOutput) {
96 delete fOutput;
97 fOutput = 0;
98 }
99 if (fVHF) {
100 delete fVHF;
101 fVHF = 0;
102 }
82bb8d4b 103
104 // if(fArrayBinLimits) {
105 //delete fArrayBinLimits;
106 //fArrayBinLimits= 0;
107 //}
108
d486095a 109}
82bb8d4b 110//_________________________________________________________________
111void AliAnalysisTaskSEDplus::SetMassLimits(Float_t range){
112 fUpmasslimit = 1.865+range;
113 fLowmasslimit = 1.865-range;
114}
115//_________________________________________________________________
116void AliAnalysisTaskSEDplus::SetMassLimits(Float_t lowlimit, Float_t uplimit){
117 if(uplimit>lowlimit)
118 {
119 fUpmasslimit = lowlimit;
120 fLowmasslimit = uplimit;
121 }
122}
123
124
125//________________________________________________________________________
126void AliAnalysisTaskSEDplus::SetPtBinLimit(Int_t n, Double_t* lim){
127 // define pt bins for analysis
128 if(n>kMaxPtBins){
129 printf("Max. number of Pt bins = %d\n",kMaxPtBins);
130 fNPtBins=kMaxPtBins;
131 fArrayBinLimits[0]=0.;
132 fArrayBinLimits[1]=2.;
133 fArrayBinLimits[2]=3.;
134 fArrayBinLimits[3]=5.;
135 for(Int_t i=4; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
136 }else{
137 fNPtBins=n-1;
138 fArrayBinLimits[0]=lim[0];
139 for(Int_t i=1; i<fNPtBins+1; i++)
140 if(lim[i]>fArrayBinLimits[i-1]){
141 fArrayBinLimits[i]=lim[i];
142 }
143 else {
144 fArrayBinLimits[i]=fArrayBinLimits[i-1];
145 }
146 for(Int_t i=fNPtBins; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
147 }
148 if(fDebug > 1){
149 printf("Number of Pt bins = %d\n",fNPtBins);
150 for(Int_t i=0; i<fNPtBins+1; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fArrayBinLimits[i],fArrayBinLimits[i+1]);
151 }
152}
153//_________________________________________________________________
154Double_t AliAnalysisTaskSEDplus::GetPtBinLimit(Int_t ibin){
155 if(ibin>fNPtBins)return -1;
156 return fArrayBinLimits[ibin];
157}
158
159//_________________________________________________________________
160void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS){
161
162/*
163 * Fill the Like Sign histograms
164 */
165
166 Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};
167
168 //count pos/neg tracks
169 Int_t nPosTrks=0,nNegTrks=0;
170 //counter for particles passing single particle cuts
171 Int_t nspcplus=0;
172 Int_t nspcminus=0;
173
174 for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
175 AliAODTrack *track = aod->GetTrack(it);
176 if(track->Charge()>0){
177 nPosTrks++;
178 if(track->Pt()>=0.4){
179 nspcplus++;
180 }
181 }
182 if(track->Charge()<0)
183 {
184 nNegTrks++;
185 if(track->Pt()>=0.4){
186 nspcminus++;
187 }
188 }
189 }
190
191 Int_t nOStriplets = arrayOppositeSign->GetEntriesFast();
192
193 Int_t nDplusLS=0;
194 Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
195 Int_t index;
196
197 for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
198 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
199 Bool_t unsetvtx=kFALSE;
200 if(!d->GetOwnPrimaryVtx()) {
201 d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
202 unsetvtx=kTRUE;
203 }
204 if(d->SelectDplus(fVHF->GetDplusCuts()))nDplusLS++;
205 if(unsetvtx) d->UnsetOwnPrimaryVtx();
206 }
207
208 Float_t wei2=0;
209 if(nLikeSign!=0)wei2 = (Float_t)nOStriplets/(Float_t)nLikeSign;
210 Float_t wei3=0;
211 if(nDplusLS!=0)wei3 = (Float_t)nDplusOS/(Float_t)nDplusLS;
212
213 // loop over like sign candidates
214 for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
215 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
216 Bool_t unsetvtx=kFALSE;
217 if(!d->GetOwnPrimaryVtx()) {
218 d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
219 unsetvtx=kTRUE;
220 }
221
222 if(d->SelectDplus(fVHF->GetDplusCuts())){
223
224 //set tight cuts values
225 Int_t iPtBin=-1;
226 Double_t ptCand = d->Pt();
227 if(ptCand<2.){
228 //iPtBin=0;
229 cutsDplus[7]=0.08;
230 cutsDplus[8]=0.5;
13808a30 231 cutsDplus[9]=0.979;
82bb8d4b 232 cutsDplus[10]=0.0055;
233 }
234 else if(ptCand>2. && ptCand<3){
235 //iPtBin=1;
236 cutsDplus[7]=0.08;
237 cutsDplus[8]=0.5;
13808a30 238 cutsDplus[9]=0.991;
82bb8d4b 239 cutsDplus[10]=0.005;
240 }else if(ptCand>3. && ptCand<5){
241 //iPtBin=2;
242 cutsDplus[7]=0.1;
243 cutsDplus[8]=0.5;
13808a30 244 cutsDplus[9]=0.995;
82bb8d4b 245 cutsDplus[10]=0.0035;
246 }else{
247 //iPtBin=3;
248 cutsDplus[7]=0.1;
249 cutsDplus[8]=0.5;
13808a30 250 cutsDplus[9]=0.997;
82bb8d4b 251 cutsDplus[10]=0.001;
252 }
253
254 for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
255 if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
256 }
257
258 if(iPtBin<0){
259 return;
260 }
261
262 Bool_t passTightCuts=d->SelectDplus(cutsDplus);
263
264 Int_t sign= d->GetCharge();
265 Float_t wei=1;
266 Float_t wei4=1;
267 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
268
269 wei=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.);
270 wei4=3.*(Float_t)nspcminus/((Float_t)nspcplus-2.);
271 }
272
273 if(sign<0&&nNegTrks>2&&nspcminus>2){
274 wei=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.);
275 wei4=3.*(Float_t)nspcplus/((Float_t)nspcminus-2.);
276
277 }
278
279 Float_t invMass = d->InvMassDplus();
280
281
282 index=GetLSHistoIndex(iPtBin);
283 fMassHistLS[index]->Fill(invMass,wei);
284 fMassHistLS[index+1]->Fill(invMass);
285 fMassHistLS[index+2]->Fill(invMass,wei2);
286 fMassHistLS[index+3]->Fill(invMass,wei3);
287 fMassHistLS[index+4]->Fill(invMass,wei4);
288
289
290 if(passTightCuts){
291 fMassHistLSTC[index]->Fill(invMass,wei);
292 fMassHistLSTC[index+1]->Fill(invMass);
293 fMassHistLSTC[index+2]->Fill(invMass,wei2);
294 fMassHistLSTC[index+3]->Fill(invMass,wei3);
295 fMassHistLSTC[index+4]->Fill(invMass,wei4);
296 }
297 }
298 if(unsetvtx) d->UnsetOwnPrimaryVtx();
299 }
300
301 //printf("------------ N. of positive tracks in Event ----- %d \n", nPosTrks);
302 //printf("------------ N. of negative tracks in Event ----- %d \n", nNegTrks);
303
304 // printf("LS analysis...done\n");
305
306}
d486095a 307
308//________________________________________________________________________
309void AliAnalysisTaskSEDplus::Init()
310{
311 // Initialization
312
313 if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
314
315 gROOT->LoadMacro("ConfigVertexingHF.C");
316
317 fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
318 fVHF->PrintStatus();
319
320 return;
321}
322
323//________________________________________________________________________
324void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
325{
326 // Create the output container
327 //
328 if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");
329
330 // Several histograms are more conveniently managed in a TList
331 fOutput = new TList();
332 fOutput->SetOwner();
1f4e9722 333 fOutput->SetName("OutputHistos");
334
1f4e9722 335 TString hisname;
82bb8d4b 336 Int_t index=0;
337 Int_t indexLS=0;
338 for(Int_t i=0;i<fNPtBins;i++){
1f4e9722 339
82bb8d4b 340 index=GetHistoIndex(i);
341 indexLS=GetLSHistoIndex(i);
342 hisname.Form("hMassPt%d",i);
343 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
344 fMassHist[index]->Sumw2();
1f4e9722 345 hisname.Form("hMassPt%dTC",i);
82bb8d4b 346 fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
347 fMassHistTC[index]->Sumw2();
348 hisname.Form("hLSPt%dTC",i);
349 fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
350 fMassHistLSTC[indexLS]->Sumw2();
351 hisname.Form("hLSPt%dLC",i);
352 fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
353 fMassHistLS[indexLS]->Sumw2();
354
355 index=GetSignalHistoIndex(i);
356 indexLS++;
357 hisname.Form("hSigPt%d",i);
358 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
359 fMassHist[index]->Sumw2();
1f4e9722 360 hisname.Form("hSigPt%dTC",i);
82bb8d4b 361 fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
362 fMassHistTC[index]->Sumw2();
363 hisname.Form("hLSPt%dLCnw",i);
364 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
365 fMassHistLS[indexLS]->Sumw2();
366 hisname.Form("hLSPt%dTCnw",i);
367 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
368 fMassHistLSTC[indexLS]->Sumw2();
369
370 index=GetBackgroundHistoIndex(i);
371 indexLS++;
372 hisname.Form("hBkgPt%d",i);
373 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
374 fMassHist[index]->Sumw2();
1f4e9722 375 hisname.Form("hBkgPt%dTC",i);
82bb8d4b 376 fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
377 fMassHistTC[index]->Sumw2();
378 hisname.Form("hLSPt%dLCntrip",i);
379 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
380 fMassHistLS[indexLS]->Sumw2();
381 hisname.Form("hLSPt%dTCntrip",i);
382 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
383 fMassHistLSTC[indexLS]->Sumw2();
384
385 indexLS++;
386 hisname.Form("hLSPt%dLCntripsinglecut",i);
387 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
388 fMassHistLS[indexLS]->Sumw2();
389 hisname.Form("hLSPt%dTCntripsinglecut",i);
390 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
391 fMassHistLSTC[indexLS]->Sumw2();
392
393 indexLS++;
394 hisname.Form("hLSPt%dLCspc",i);
395 fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
396 fMassHistLS[indexLS]->Sumw2();
397 hisname.Form("hLSPt%dTCspc",i);
398 fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
399 fMassHistLSTC[indexLS]->Sumw2();
1f4e9722 400 }
8c34bd86 401
82bb8d4b 402 for(Int_t i=0; i<3*fNPtBins; i++){
403 fOutput->Add(fMassHist[i]);
404 fOutput->Add(fMassHistTC[i]);
405 }
406 for(Int_t i=0; i<5*fNPtBins&&fDoLS; i++){
407 fOutput->Add(fMassHistLS[i]);
408 fOutput->Add(fMassHistLSTC[i]);
409 }
4afc48a2 410
4afc48a2 411
82bb8d4b 412 fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
413 fHistNEvents->Sumw2();
414 fHistNEvents->SetMinimum(0);
415 fOutput->Add(fHistNEvents);
4afc48a2 416
1f4e9722 417 if(fFillNtuple){
418 OpenFile(2); // 2 is the slot number of the ntuple
4afc48a2 419 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");
1f4e9722 420 }
421
d486095a 422 return;
423}
424
425//________________________________________________________________________
426void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
427{
428 // Execute analysis for current event:
429 // heavy flavor candidates association to MC truth
430
82bb8d4b 431 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
432 fHistNEvents->Fill(0); // count event
4afc48a2 433 // Post the data already here
434 PostData(1,fOutput);
435
b557eb43 436 TClonesArray *array3Prong = 0;
4afc48a2 437 TClonesArray *arrayLikeSign =0;
b557eb43 438 if(!aod && AODEvent() && IsStandardAOD()) {
439 // In case there is an AOD handler writing a standard AOD, use the AOD
440 // event in memory rather than the input (ESD) event.
441 aod = dynamic_cast<AliAODEvent*> (AODEvent());
442 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
443 // have to taken from the AOD event hold by the AliAODExtension
444 AliAODHandler* aodHandler = (AliAODHandler*)
445 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
446 if(aodHandler->GetExtensions()) {
447 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
448 AliAODEvent *aodFromExt = ext->GetAOD();
449 array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
4afc48a2 450 arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
b557eb43 451 }
452 } else {
453 array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
4afc48a2 454 arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
b557eb43 455 }
8931c313 456
d486095a 457 if(!array3Prong) {
458 printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
459 return;
460 }
4afc48a2 461 if(!arrayLikeSign) {
462 printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
82bb8d4b 463 return;
4afc48a2 464 }
465
466
467 TClonesArray *arrayMC=0;
468 AliAODMCHeader *mcHeader=0;
d486095a 469
470 // AOD primary vertex
1f4e9722 471 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
472 // vtx1->Print();
473
d486095a 474 // load MC particles
4afc48a2 475 if(fReadMC){
476
477 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
478 if(!arrayMC) {
479 printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
82bb8d4b 480 // return;
4afc48a2 481 }
1f4e9722 482
d486095a 483 // load MC header
4afc48a2 484 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
485 if(!mcHeader) {
82bb8d4b 486 printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
487 return;
4afc48a2 488 }
d486095a 489 }
4afc48a2 490
1f4e9722 491 Int_t n3Prong = array3Prong->GetEntriesFast();
492 if(fDebug>1) printf("Number of D+->Kpipi: %d\n",n3Prong);
4afc48a2 493
494
495 Int_t nOS=0;
82bb8d4b 496 Int_t index;
1f4e9722 497 Int_t pdgDgDplustoKpipi[3]={321,211,211};
498 Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};
499 for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
500 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
501
502
503 Bool_t unsetvtx=kFALSE;
504 if(!d->GetOwnPrimaryVtx()){
505 d->SetOwnPrimaryVtx(vtx1);
506 unsetvtx=kTRUE;
507 }
4afc48a2 508
1f4e9722 509 if(d->SelectDplus(fVHF->GetDplusCuts())) {
82bb8d4b 510 Int_t iPtBin = -1;
1f4e9722 511 Double_t ptCand = d->Pt();
82bb8d4b 512
1f4e9722 513 if(ptCand<2.){
82bb8d4b 514 //iPtBin=0;
1f4e9722 515 cutsDplus[7]=0.08;
516 cutsDplus[8]=0.5;
13808a30 517 cutsDplus[9]=0.979;
82bb8d4b 518 cutsDplus[10]=0.0055;
1f4e9722 519 }
520 else if(ptCand>2. && ptCand<3){
82bb8d4b 521 //iPtBin=1;
1f4e9722 522 cutsDplus[7]=0.08;
523 cutsDplus[8]=0.5;
13808a30 524 cutsDplus[9]=0.991;
82bb8d4b 525 cutsDplus[10]=0.005;
1f4e9722 526 }else if(ptCand>3. && ptCand<5){
82bb8d4b 527 //iPtBin=2;
1f4e9722 528 cutsDplus[7]=0.1;
529 cutsDplus[8]=0.5;
13808a30 530 cutsDplus[9]=0.995;
82bb8d4b 531 cutsDplus[10]=0.0035;
1f4e9722 532 }else{
82bb8d4b 533 //iPtBin=3;
1f4e9722 534 cutsDplus[7]=0.1;
535 cutsDplus[8]=0.5;
13808a30 536 cutsDplus[9]=0.997;
82bb8d4b 537 cutsDplus[10]=0.001;
538 }
539
540 for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
541 if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
1f4e9722 542 }
82bb8d4b 543
1f4e9722 544 Bool_t passTightCuts=d->SelectDplus(cutsDplus);
4afc48a2 545 Int_t labDp=-1;
546 Float_t deltaPx=0.;
547 Float_t deltaPy=0.;
548 Float_t deltaPz=0.;
549 Float_t truePt=0.;
550 Float_t xDecay=0.;
551 Float_t yDecay=0.;
552 Float_t zDecay=0.;
553 Float_t pdgCode=-2;
554 if(fReadMC){
555 labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
556 if(labDp>=0){
557 AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
558 AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
559 deltaPx=partDp->Px()-d->Px();
560 deltaPy=partDp->Py()-d->Py();
561 deltaPz=partDp->Pz()-d->Pz();
562 truePt=partDp->Pt();
563 xDecay=dg0->Xv();
564 yDecay=dg0->Yv();
565 zDecay=dg0->Zv();
566 pdgCode=TMath::Abs(partDp->GetPdgCode());
567 }else{
568 pdgCode=-1;
569 }
570 }
1f4e9722 571 Double_t invMass=d->InvMassDplus();
572
4afc48a2 573 Float_t tmp[22];
574 if(fFillNtuple){
575 tmp[0]=pdgCode;
576 tmp[1]=deltaPx;
577 tmp[2]=deltaPy;
578 tmp[3]=deltaPz;
579 tmp[4]=truePt;
580 tmp[5]=xDecay;
581 tmp[6]=yDecay;
582 tmp[7]=zDecay;
583 tmp[8]=d->PtProng(0);
584 tmp[9]=d->PtProng(1);
585 tmp[10]=d->PtProng(2);
586 tmp[11]=d->Pt();
587 tmp[12]=d->CosPointingAngle();
588 tmp[13]=d->DecayLength();
589 tmp[14]=d->Xv();
590 tmp[15]=d->Yv();
591 tmp[16]=d->Zv();
592 tmp[17]=d->InvMassDplus();
593 tmp[18]=d->GetSigmaVert();
594 tmp[19]=d->Getd0Prong(0);
595 tmp[20]=d->Getd0Prong(1);
596 tmp[21]=d->Getd0Prong(2);
597 fNtupleDplus->Fill(tmp);
598 PostData(2,fNtupleDplus);
599 }
d486095a 600
82bb8d4b 601 if(iPtBin>=0){
602
603 index=GetHistoIndex(iPtBin);
604 fMassHist[index]->Fill(invMass);
605 if(passTightCuts){
606 fMassHistTC[index]->Fill(invMass);
607
608 }
609
610 if(fReadMC){
611 if(labDp>=0) {
612 index=GetSignalHistoIndex(iPtBin);
613 fMassHist[index]->Fill(invMass);
614
615 if(passTightCuts){
616 fMassHistTC[index]->Fill(invMass);
617
618 }
619
620 }else{
621 index=GetBackgroundHistoIndex(iPtBin);
622 fMassHist[index]->Fill(invMass);
623
624 if(passTightCuts){
625 fMassHistTC[index]->Fill(invMass);
626
627 }
4afc48a2 628 }
fc8d975b 629 }
8c34bd86 630 }
82bb8d4b 631 /*
632 //start OS analysis
633 if(labDp<0)fHistOSbkg->Fill(d->InvMassDplus());
634 fHistOS->Fill(d->InvMassDplus());
635 */
4afc48a2 636 nOS++;
d486095a 637 }
1f4e9722 638 if(unsetvtx) d->UnsetOwnPrimaryVtx();
639 }
4afc48a2 640
641 //start LS analysis
642 if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
82bb8d4b 643
4afc48a2 644 PostData(1,fOutput);
d486095a 645 return;
646}
647
4afc48a2 648
649
82bb8d4b 650//________________________________________________________________________
d486095a 651void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
652{
653 // Terminate analysis
654 //
655 if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
82bb8d4b 656
d486095a 657 fOutput = dynamic_cast<TList*> (GetOutputData(1));
658 if (!fOutput) {
659 printf("ERROR: fOutput not available\n");
660 return;
661 }
82bb8d4b 662 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
663
664 TString hisname;
665 Int_t index=0;
666 Int_t indexLS=0;
667 for(Int_t i=0;i<fNPtBins;i++){
668 index=GetHistoIndex(i);
669 if(fDoLS)indexLS=GetLSHistoIndex(i);
670 hisname.Form("hMassPt%d",i);
671 fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
672
673 hisname.Form("hMassPt%dTC",i);
674 fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
675 if(fDoLS){
676 hisname.Form("hLSPt%dTC",i);
677 fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
678
679 hisname.Form("hLSPt%dLC",i);
680 fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
681 }
682
683 index=GetSignalHistoIndex(i);
684 if(fDoLS)indexLS++;
685 hisname.Form("hSigPt%d",i);
686 fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
687
688 hisname.Form("hSigPt%dTC",i);
689 fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
690 if(fDoLS){
691 hisname.Form("hLSPt%dLCnw",i);
692 fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
693
694 hisname.Form("hLSPt%dTCnw",i);
695 fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
696 }
697
698 index=GetBackgroundHistoIndex(i);
699 if(fDoLS)indexLS++;
700 hisname.Form("hBkgPt%d",i);
701 fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
702
703 hisname.Form("hBkgPt%dTC",i);
704 fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
705 if(fDoLS){
706 hisname.Form("hLSPt%dLCntrip",i);
707 fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
708
709 hisname.Form("hLSPt%dTCntrip",i);
710 fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
711
712
713 indexLS++;
714 hisname.Form("hLSPt%dLCntripsinglecut",i);
715 fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
716
717 hisname.Form("hLSPt%dTCntripsinglecut",i);
718 fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
719
720
721 indexLS++;
722 hisname.Form("hLSPt%dLCspc",i);
723 fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
724
725 hisname.Form("hLSPt%dTCspc",i);
726 fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
727 }
728
729 }
730
1f4e9722 731 if(fFillNtuple){
732 fNtupleDplus = dynamic_cast<TNtuple*>(GetOutputData(2));
4afc48a2 733 }
fc8d975b 734
13808a30 735 TCanvas *c1=new TCanvas("c1","D+ invariant mass distribution",500,500);
736 c1->cd();
737 TH1F *hMassPt=(TH1F*)fOutput->FindObject("hMassPt3TC");
738 hMassPt->SetLineColor(kBlue);
739 hMassPt->SetXTitle("M[GeV/c^{2}]");
740 hMassPt->Draw();
741
742 return;
d486095a 743}