]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliAnalysisTaskSEDStarSpectra.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New class for PID of HF candidates (R. Romita)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDStarSpectra.cxx
CommitLineData
645e004b 1/**************************************************************************
2 * Copyright(c) 1998-2009, 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// Base class for DStar Analysis
18//
19// Side Band and like sign background are implemented in the macro
20//
21// Utrech Optimized cuts used and PID is on request (flag in che .C)
22// The D* spectra study is done in pt bins:
23//
24// [0,1] [1,2] [2,3] [3,5] [5,8] [8,14]
25//
26//-----------------------------------------------------------------------
27//
28// Author A.Grelli
29// ERC-QGP Utrecht University - a.grelli@uu.nl
30// Contributors: C.Ivan - Utrecht University
31//
32//-----------------------------------------------------------------------
33
34#include <TSystem.h>
35#include <TParticle.h>
36#include <TH1I.h>
37#include "TROOT.h"
38
39#include "AliPID.h"
40#include "AliTPCPIDResponse.h"
41#include "AliStack.h"
42#include "AliMCEvent.h"
43#include "AliAnalysisManager.h"
44#include "AliAODMCHeader.h"
45#include "AliAODHandler.h"
46#include "AliLog.h"
47#include "AliAODVertex.h"
48#include "AliAODJet.h"
49#include "AliAODRecoDecay.h"
50#include "AliAODRecoDecayHF.h"
51#include "AliAODRecoCascadeHF.h"
52#include "AliAODRecoDecayHF2Prong.h"
53#include "AliAnalysisVertexingHF.h"
54#include "AliESDtrack.h"
55#include "AliAODMCParticle.h"
56#include "AliAnalysisTaskSEDStarSpectra.h"
57
58ClassImp(AliAnalysisTaskSEDStarSpectra)
59
60//__________________________________________________________________________
61AliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra():
62 AliAnalysisTaskSE(),
63 fEvents(0),
64 fVHF(0),
65 fMinITSClusters(0),
66 fMinITSClustersSoft(0),
67 fUseMCInfo(kTRUE),
68 fOutput(0),
69 fNSigma(3),
70 fPID(kTRUE),
71 fAODTrack(0),
72 fMCDStarPt(0),
73 fCEvents(0),
74 fDStarMass(0),
75 fTrueDiff(0),
76 fTrueDiff2(0),
77 fInvMass(0),
78 fInvMass1(0),
79 fInvMass2(0),
80 fInvMass3(0),
81 fInvMass4(0),
82 fInvMass5(0),
83 fPtDStar(0),
84 fDStar(0),
85 fDiff(0),
86 fDiff1(0),
87 fDiff2(0),
88 fDiff3(0),
89 fDiff4(0),
90 fDiff5(0),
91 fDiffSideBand(0),
92 fDiffSideBand1(0),
93 fDiffSideBand2(0),
94 fDiffSideBand3(0),
95 fDiffSideBand4(0),
96 fDiffSideBand5(0),
97 fDiffWrongSign(0),
98 fDiffWrongSign1(0),
99 fDiffWrongSign2(0),
100 fDiffWrongSign3(0),
101 fDiffWrongSign4(0),
102 fDiffWrongSign5(0)
103
104{
105 //
106 // Default ctor
107 //
108}
109//___________________________________________________________________________
110AliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra(const Char_t* name) :
111 AliAnalysisTaskSE(name),
112 fEvents(0),
113 fVHF(0),
114 fMinITSClusters(0),
115 fMinITSClustersSoft(0),
116 fUseMCInfo(kTRUE),
117 fOutput(0),
118 fNSigma(3),
119 fPID(kTRUE),
120 fAODTrack(0),
121 fMCDStarPt(0),
122 fCEvents(0),
123 fDStarMass(0),
124 fTrueDiff(0),
125 fTrueDiff2(0),
126 fInvMass(0),
127 fInvMass1(0),
128 fInvMass2(0),
129 fInvMass3(0),
130 fInvMass4(0),
131 fInvMass5(0),
132 fPtDStar(0),
133 fDStar(0),
134 fDiff(0),
135 fDiff1(0),
136 fDiff2(0),
137 fDiff3(0),
138 fDiff4(0),
139 fDiff5(0),
140 fDiffSideBand(0),
141 fDiffSideBand1(0),
142 fDiffSideBand2(0),
143 fDiffSideBand3(0),
144 fDiffSideBand4(0),
145 fDiffSideBand5(0),
146 fDiffWrongSign(0),
147 fDiffWrongSign1(0),
148 fDiffWrongSign2(0),
149 fDiffWrongSign3(0),
150 fDiffWrongSign4(0),
151 fDiffWrongSign5(0)
152{
153 //
154 // Constructor. Initialization of Inputs and Outputs
155 //
156 Info("AliAnalysisTaskSEDStarSpectra","Calling Constructor");
157
158 DefineOutput(1,TList::Class());
159}
160
161//___________________________________________________________________________
162AliAnalysisTaskSEDStarSpectra& AliAnalysisTaskSEDStarSpectra::operator=(const AliAnalysisTaskSEDStarSpectra& c)
163{
164 //
165 // Assignment operator
166 //
167 if (this!=&c) {
168 AliAnalysisTaskSE::operator=(c) ;
169 }
170 return *this;
171}
172
173//___________________________________________________________________________
174AliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra(const AliAnalysisTaskSEDStarSpectra& c) :
175 AliAnalysisTaskSE(c),
176 fEvents(c.fEvents),
177 fVHF(c.fVHF),
178 fMinITSClusters(c.fMinITSClusters),
179 fMinITSClustersSoft(c.fMinITSClustersSoft),
180 fUseMCInfo(c.fUseMCInfo),
181 fOutput(c.fOutput),
182 fNSigma(c.fNSigma),
183 fPID(c.fPID),
184 fAODTrack(c.fAODTrack),
185 fMCDStarPt(c.fMCDStarPt),
186 fCEvents(c.fCEvents),
187 fDStarMass(c.fDStarMass),
188 fTrueDiff(c.fTrueDiff),
189 fTrueDiff2(c.fTrueDiff2),
190 fInvMass(c.fInvMass),
191 fInvMass1(c.fInvMass1),
192 fInvMass2(c.fInvMass2),
193 fInvMass3(c.fInvMass3),
194 fInvMass4(c.fInvMass4),
195 fInvMass5(c.fInvMass5),
196 fPtDStar(c.fPtDStar),
197 fDStar(c.fDStar),
198 fDiff(c.fDiff),
199 fDiff1(c.fDiff1),
200 fDiff2(c.fDiff2),
201 fDiff3(c.fDiff3),
202 fDiff4(c.fDiff4),
203 fDiff5(c.fDiff5),
204 fDiffSideBand(c.fDiffSideBand),
205 fDiffSideBand1(c.fDiffSideBand1),
206 fDiffSideBand2(c.fDiffSideBand2),
207 fDiffSideBand3(c.fDiffSideBand3),
208 fDiffSideBand4(c.fDiffSideBand4),
209 fDiffSideBand5(c.fDiffSideBand5),
210 fDiffWrongSign(c.fDiffWrongSign),
211 fDiffWrongSign1(c.fDiffWrongSign1),
212 fDiffWrongSign2(c.fDiffWrongSign2),
213 fDiffWrongSign3(c.fDiffWrongSign3),
214 fDiffWrongSign4(c.fDiffWrongSign4),
215 fDiffWrongSign5(c.fDiffWrongSign5)
216{
217 //
218 // Copy Constructor
219 //
220}
221
222//___________________________________________________________________________
223AliAnalysisTaskSEDStarSpectra::~AliAnalysisTaskSEDStarSpectra() {
224 //
225 // destructor
226 //
227 Info("~AliAnalysisTaskSEDStarSpectra","Calling Destructor");
228
229 if (fOutput) {
230 delete fOutput;
231 fOutput = 0;
232 }
233 if (fVHF) {
234 delete fVHF;
235 fVHF = 0;
236 }
237}
238//_________________________________________________
239void AliAnalysisTaskSEDStarSpectra::Init(){
240 //
241 // Initialization
242 //
243
244 if(fDebug > 1) printf("AnalysisTaskSEDStarSpectra::Init() \n");
245
246 gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/ConfigVertexingHF.C");
247 fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
248 //fVHF->PrintStatus();
249
250 return;
251}
252
253//_________________________________________________
254void AliAnalysisTaskSEDStarSpectra::UserExec(Option_t *)
255{
256 // user exec
257 if (!fInputEvent) {
258 Error("UserExec","NO EVENT FOUND!");
259 return;
260 }
261
262 fCEvents->Fill(1);
263 // Load the event
264 fEvents++;
265 AliInfo(Form("Event %d",fEvents));
266 if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
267 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
268 TClonesArray *arrayDStartoD0pi=0;
269 Init();
270 if(!aodEvent && AODEvent() && IsStandardAOD()) {
271 // In case there is an AOD handler writing a standard AOD, use the AOD
272 // event in memory rather than the input (ESD) event.
273 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
274 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
275 // have to taken from the AOD event hold by the AliAODExtension
276 AliAODHandler* aodHandler = (AliAODHandler*)
277 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
278 if(aodHandler->GetExtensions()) {
279 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
280 AliAODEvent *aodFromExt = ext->GetAOD();
281 arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
282 }
283 } else {
284 arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
285 }
286
287 // AOD primary vertex
288 AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
289
290 // counters for efficiencies
291 Int_t icountReco = 0;
292
293 //D* and D0 prongs needed to MatchToMC method
294 Int_t pdgDgDStartoD0pi[2]={421,211};
295 Int_t pdgDgD0toKpi[2]={321,211};
296
297 if (!arrayDStartoD0pi){
298 AliInfo("Could not find array of HF vertices, skipping the event");
299 return;
300 }else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
301
302 // loop over the tracks to search for candidates soft pion
303
304 for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
305
306 // D* candidates
307 AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
308
309 // D0 from the reco cascade
310 AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
311 Bool_t unsetvtx=kFALSE;
312
313 Double_t finvM =0;
314 Double_t finvMDStar = 0;
315 // needed for pointing angle
316 if(!theD0particle->GetOwnPrimaryVtx()) {
317 theD0particle->SetOwnPrimaryVtx(vtx1);
318 unsetvtx=kTRUE;
319 }
320 Bool_t isDStar = 0;
321
322 // mc analysis
323 if(fUseMCInfo){
324 //MC array need for maching
325 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
326 if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
327 // find associated MC particle for D* ->D0toKpi
328 Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
329 if(mcLabel>=0) isDStar = 1;
330 }
331
332 // soft pion
333 AliAODTrack *track2 = (AliAODTrack*)dstarD0pi->GetBachelor();
334
335 //D0tokpi
336 AliAODTrack *track0 = (AliAODTrack*)theD0particle->GetDaughter(0);
337 AliAODTrack *track1 = (AliAODTrack*)theD0particle->GetDaughter(1);
338
339 Double_t pt = dstarD0pi->Pt();
340
341 //kaon PID on low pt D*
342 if(pt<3 && fPID){
343 if (dstarD0pi->Charge()>0){
344 if(!SelectPID(track1,fNSigma)) continue;
345 }else{
346 if(!SelectPID(track0,fNSigma)) continue;
347 }
348 }
349
350 // reft in ITS for soft pion
351 //if((!(track2->GetStatus()&AliESDtrack::kITSrefit))) continue;
352
353 // cut in acceptance for the soft pion and for the D0 daughters
354 Bool_t acceptanceProng0 = (TMath::Abs(theD0particle->EtaProng(0))<= 0.9 && theD0particle->PtProng(0) >= 0.1);
355 Bool_t acceptanceProng1 = (TMath::Abs(theD0particle->EtaProng(1))<= 0.9 && theD0particle->PtProng(1) >= 0.1);
356 // soft pion acceptance ... is it fine 0.9?????
357 Bool_t acceptanceProng2 = (TMath::Abs(track2->Eta())<= 1.0 && track2->Pt() >= 0.05);
358
359 if (acceptanceProng0 && acceptanceProng1 && acceptanceProng2) {
360 AliDebug(2,"D* reco daughters in acceptance");
361
362 // cut on the min n. of clusters in ITS for the D0 and soft pion
363 Int_t ncls0=0,ncls1=0,ncls2=0;
364 for(Int_t l=0;l<6;l++) {
365 if(TESTBIT(track0->GetITSClusterMap(),l)) ncls0++;
366 if(TESTBIT(track1->GetITSClusterMap(),l)) ncls1++;
367 if(TESTBIT(track2->GetITSClusterMap(),l)) ncls2++;
368 }
369 // see AddTask for soft pion and D0 prongs ITS clusters request
370 if (ncls0 >= fMinITSClusters && ncls1 >= fMinITSClusters && ncls2>=fMinITSClustersSoft) {
371
372 // tag the D* decay do not double count background
373 Int_t decayTag = track2->Charge();
374 // D0 pt needed for the cuts
375 Double_t ptD0 = theD0particle->Pt();
376
377 Int_t okD0 = 0;
378 Int_t okD0bar = 0;
379 Int_t okD0WrongSign =0;
380 Int_t okD0barWrongSign =0;
381
382 // optimized D* cuts from UU
fca7f2a2 383 Bool_t tCutDStar = SetUtrechtSelections(ptD0);
645e004b 384 Bool_t tCutOk = kFALSE;
385 if(tCutDStar) tCutOk = theD0particle->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar);
386
387 // flags for wrong sign
388 okD0WrongSign=okD0;
389 okD0barWrongSign=okD0bar;
390
391 // search for D*
392 if(tCutOk){
393 // correct decay
fca7f2a2 394 if(decayTag>0) okD0bar = 0;
395 if(decayTag<0) okD0 = 0;
645e004b 396
397 finvM = dstarD0pi->InvMassD0();
398
399 if(okD0 == 1) fInvMass->Fill(finvM);
400 if(okD0bar == 1) fInvMass->Fill(finvM);
401
402 //DStar invariant mass
403 finvMDStar = dstarD0pi->InvMassDstarKpipi();
404
405 if(finvM >= 1.829 && finvM <= 1.901){ // ~3 sigma cut on D0 mass
406
407 if(isDStar == 1) {
408 fDStarMass->Fill(finvMDStar);
409 fTrueDiff ->Fill(dstarD0pi->DeltaInvMass());
410 fTrueDiff2->Fill(pt,dstarD0pi->DeltaInvMass());
411 fMCDStarPt->Fill(pt);
412 }
413
414 // D* candidates - pt integrated
415 if(okD0==1){
416 fDStar->Fill(finvMDStar);
417 fDiff->Fill(dstarD0pi->DeltaInvMass()); // M(Kpipi)-M(Kpi)
418 }else if(okD0bar==1){
419 fDStar->Fill(finvMDStar);
420 fDiff->Fill(dstarD0pi->DeltaInvMass()); // M(Kpipi)-M(Kpi)
421 }
422
423 //D* candidates - pt bins
424 if(pt>1 && pt<=2){ // [1-2]
425 if(okD0==1){
426 fDiff1->Fill(dstarD0pi->DeltaInvMass());
427 fInvMass1->Fill(finvM);
428 }else if(okD0bar==1){
429 fDiff1->Fill(dstarD0pi->DeltaInvMass());
430 fInvMass1->Fill(finvM);
431 }
432 }
433 if( pt>2 && pt<=3){ // [2-3]
434 if(okD0==1){
435 fDiff2->Fill(dstarD0pi->DeltaInvMass());
436 fInvMass2->Fill(finvM);
437 }else if(okD0bar==1){
438 fDiff2->Fill(dstarD0pi->DeltaInvMass());
439 fInvMass2->Fill(finvM);
440 }
441 }
442 if(pt>3 && pt<=5){ // [3-5]
443 if(okD0==1){
444 fDiff3->Fill(dstarD0pi->DeltaInvMass());
445 fInvMass3->Fill(finvM);
446 }else if(okD0bar==1){
447 fDiff3->Fill(dstarD0pi->DeltaInvMass());
448 fInvMass3->Fill(finvM);
449 }
450 }
451 if( pt>5 && pt<8){ // [5-8]
452 if(okD0==1){
453 fDiff4->Fill(dstarD0pi->DeltaInvMass());
454 fInvMass4->Fill(finvM);
455 }else if(okD0bar==1){
456 fDiff4->Fill(dstarD0pi->DeltaInvMass());
457 fInvMass4->Fill(finvM);
458 }
459 }
460 if(pt>=8){ // [>8]
461 if(okD0==1){
462 fDiff5->Fill(dstarD0pi->DeltaInvMass());
463 fInvMass5->Fill(finvM);
464 }else if(okD0bar==1){
465 fDiff5->Fill(dstarD0pi->DeltaInvMass());
466 fInvMass5->Fill(finvM);
467 }
468 }
469
470 // D* pt distribution
471 if((dstarD0pi->DeltaInvMass())>=0.143920 && (dstarD0pi->DeltaInvMass())<=0.14702){
472 if(okD0==1) fPtDStar->Fill(pt);
473 if(okD0bar==1) fPtDStar->Fill(pt);
474 }
475
476 }
477 // evaluate side band background
478 SideBandBackground(finvM,finvMDStar,pt,okD0,okD0bar);
479
480 finvM = 0;
481 finvMDStar = 0;
482
483 } // end cuts
484
485 // wrong sign background
486 if(tCutOk){
487 // correct decay
fca7f2a2 488 if(decayTag>0) okD0WrongSign = 0;
489 if(decayTag<0) okD0barWrongSign = 0;
490
645e004b 491 //wrong D0 inv mass
492 if(okD0WrongSign==1){
493 finvM = theD0particle->InvMassD0();
494 }else if(okD0WrongSign==0){
495 finvM = theD0particle->InvMassD0bar();
496 }
497 // wrong D* inv mass
498 Double_t px[3],py[3],pz[3];
499 UInt_t pdg[3]={321,211,211};
500 pdg[0] = (decayTag>0 ? 321 : 211); // positive daughter of D0
501 px[0] = theD0particle->PxProng(0);
502 py[0] = theD0particle->PyProng(0);
503 pz[0] = theD0particle->PzProng(0);
504 pdg[1] = (decayTag>0 ? 211 : 321); // negative daughter of D0
505 px[1] = theD0particle->PxProng(1);
506 py[1] = theD0particle->PyProng(1);
507 pz[1] = theD0particle->PzProng(1);
508 pdg[2] = 211; // soft pion
509 px[2] = track2->Px();
510 py[2] = track2->Py();
511 pz[2] = track2->Pz();
512
513 Short_t dummycharge=0;
514 Double_t dummyd0[3]={0,0,0};
515 AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,3,dummycharge,px,py,pz,dummyd0);
516
517 finvMDStar = rd->InvMass(3,pdg);
518
519 delete rd; rd=NULL;
520
521 if(finvM >= 1.829 && finvM <= 1.901){ // ~3 sigma cut on D0 mass
522 WrongSignForDStar(finvM,finvMDStar,pt,okD0WrongSign,okD0barWrongSign);
523 }
524 }// end of wrong sign
525 }
526 }
527 }
528
529 AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
530
531 PostData(1,fOutput);
532}
533//________________________________________ terminate ___________________________
534void AliAnalysisTaskSEDStarSpectra::Terminate(Option_t*)
535{
536 // The Terminate() function is the last function to be called during
537 // a query. It always runs on the client, it can be used to present
538 // the results graphically or save the results to file.
539
540 Info("Terminate","");
541 AliAnalysisTaskSE::Terminate();
542
543 fOutput = dynamic_cast<TList*> (GetOutputData(1));
544 if (!fOutput) {
545 printf("ERROR: fOutput not available\n");
546 return;
547 }
548
549 fMCDStarPt = dynamic_cast<TH1F*>(fOutput->FindObject("fMCDStarPt"));
550 fCEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fCEvents"));
551 fDStarMass = dynamic_cast<TH1F*>(fOutput->FindObject("fDStarMass"));
552 fTrueDiff = dynamic_cast<TH1F*>(fOutput->FindObject("fTrueDiff"));
553 fTrueDiff2 = dynamic_cast<TH2F*>(fOutput->FindObject("fTrueDiff2"));
554 fInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass"));
555 fInvMass1 = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass1"));
556 fInvMass2 = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass2"));
557 fInvMass3 = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass3"));
558 fInvMass4 = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass4"));
559 fInvMass5 = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass5"));
560 fPtDStar = dynamic_cast<TH1F*>(fOutput->FindObject("fPtDStar "));
561 fDStar = dynamic_cast<TH1F*>(fOutput->FindObject("fDStar"));
562 fDiff = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff"));
563 fDiff1 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff1"));
564 fDiff2 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff2"));
565 fDiff3 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff3"));
566 fDiff4 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff4"));
567 fDiff5 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff5"));
568 fDiffSideBand = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand"));
569 fDiffSideBand1 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand1"));
570 fDiffSideBand2 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand2"));
571 fDiffSideBand3 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand3"));
572 fDiffSideBand4 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand4"));
573 fDiffSideBand5 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand5"));
574 fDiffWrongSign = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffWrongSign"));
575 fDiffWrongSign1 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffWrongSign1"));
576 fDiffWrongSign2 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffWrongSign2"));
577 fDiffWrongSign3 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffWrongSign3"));
578 fDiffWrongSign4 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffWrongSign4"));
579 fDiffWrongSign5 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffWrongSign5"));
580
581}
582//___________________________________________________________________________
583void AliAnalysisTaskSEDStarSpectra::UserCreateOutputObjects() {
584 // output
585 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
586
587 //slot #1
588 OpenFile(1);
589 fOutput = new TList();
590 fOutput->SetOwner();
591 // define histograms
592 DefineHistoFroAnalysis();
593
594 return;
595}
596
597//___________________________________ hiostograms _______________________________________
598Bool_t AliAnalysisTaskSEDStarSpectra::DefineHistoFroAnalysis(){
599
600 // Invariant mass related histograms
601 fInvMass = new TH1F("invMass","K#pi invariant mass distribution",1500,.5,3.5);
602 fInvMass->SetStats(kTRUE);
603 fInvMass->GetXaxis()->SetTitle("M(K#pi) GeV/c");
604 fInvMass->GetYaxis()->SetTitle("Entries");
605
606 fOutput->Add(fInvMass);
607
608 fInvMass1 = new TH1F("invMass1","K#pi invariant mass distribution",3500,1.5,2.2);
609 fInvMass1->SetStats(kTRUE);
610 fInvMass1->GetXaxis()->SetTitle("M(K#pi) GeV/c");
611 fInvMass1->GetYaxis()->SetTitle("Entries");
612
613 fOutput->Add(fInvMass1);
614
615 fInvMass2 = new TH1F("invMass2","K#pi invariant mass distribution",350,1.5,2.2);
616 fInvMass2->SetStats(kTRUE);
617 fInvMass2->GetXaxis()->SetTitle("M(K#pi) GeV/c");
618 fInvMass2->GetYaxis()->SetTitle("Entries");
619
620 fOutput->Add(fInvMass2);
621
622 fInvMass3 = new TH1F("invMass3","K#pi invariant mass distribution",350,1.5,2.2);
623 fInvMass3->SetStats(kTRUE);
624 fInvMass3->GetXaxis()->SetTitle("M(K#pi) GeV/c");
625 fInvMass3->GetYaxis()->SetTitle("Entries");
626
627 fOutput->Add(fInvMass3);
628
629 fInvMass4 = new TH1F("invMass4","K#pi invariant mass distribution",350,1.5,2.2);
630 fInvMass4->SetStats(kTRUE);
631 fInvMass4->GetXaxis()->SetTitle("M(K#pi) GeV/c");
632 fInvMass4->GetYaxis()->SetTitle("Entries");
633
634 fOutput->Add(fInvMass4);
635
636 fInvMass5 = new TH1F("invMass5","K#pi invariant mass distribution",350,1.5,2.2);
637 fInvMass5->SetStats(kTRUE);
638 fInvMass5->GetXaxis()->SetTitle("M(K#pi) GeV/c");
639 fInvMass5->GetYaxis()->SetTitle("Entries");
640
641 fOutput->Add(fInvMass5);
642
643 fDStar = new TH1F("invMassDStar","DStar invariant mass after D0 cuts ",600,1.8,2.4);
644 fDStar->SetStats(kTRUE);
645 fDStar->GetXaxis()->SetTitle("GeV/c");
646 fDStar->GetYaxis()->SetTitle("Entries");
647
648 fOutput->Add(fDStar);
649
650 fCEvents = new TH1F("fCEvents","conter",10,0,10);
651 fCEvents->SetStats(kTRUE);
652 fCEvents->GetXaxis()->SetTitle("1");
653 fCEvents->GetYaxis()->SetTitle("counts");
654
655 fOutput->Add(fCEvents);
656
657 fDiff = new TH1F("Diff","M(K#pi#pi)-M(K#pi)",750,0.1,0.2);
658 fDiff->SetStats(kTRUE);
659 fDiff->GetXaxis()->SetTitle("M(K#pi#pi)-M(K#pi) GeV/c^2");
660 fDiff->GetYaxis()->SetTitle("Entries");
661
662 fOutput->Add(fDiff);
663
664 fDiff1 = new TH1F("Diff1","M(K#pi#pi)-M(K#pi)",750,0.1,0.2);
665 fDiff1->SetStats(kTRUE);
666 fDiff1->GetXaxis()->SetTitle("M(K#pi#pi)-M(K#pi) GeV/c^2");
667 fDiff1->GetYaxis()->SetTitle("Entries");
668
669 fOutput->Add(fDiff1);
670
671 fDiff2 = new TH1F("Diff2","M(K#pi#pi)-M(K#pi)",750,0.1,0.2);
672 fDiff2->SetStats(kTRUE);
673 fDiff2->GetXaxis()->SetTitle("M(K#pi#pi)-M(K#pi) GeV/c^2");
674 fDiff2->GetYaxis()->SetTitle("Entries");
675
676 fOutput->Add(fDiff2);
677
678 fDiff3 = new TH1F("Diff3","M(K#pi#pi)-M(K#pi)",750,0.1,0.2);
679 fDiff3->SetStats(kTRUE);
680 fDiff3->GetXaxis()->SetTitle("M(K#pi#pi)-M(K#pi) GeV/c^2");
681 fDiff3->GetYaxis()->SetTitle("Entries");
682
683 fOutput->Add(fDiff3);
684
685 fDiff4 = new TH1F("Diff4","M(K#pi#pi)-M(K#pi)",750,0.1,0.2);
686 fDiff4->SetStats(kTRUE);
687 fDiff4->GetXaxis()->SetTitle("M(K#pi#pi)-M(K#pi) GeV/c^2");
688 fDiff4->GetYaxis()->SetTitle("Entries");
689
690 fOutput->Add(fDiff4);
691
692 fDiff5 = new TH1F("Diff5","M(K#pi#pi)-M(K#pi)",750,0.1,0.2);
693 fDiff5->SetStats(kTRUE);
694 fDiff5->GetXaxis()->SetTitle("M(K#pi#pi)-M(K#pi) GeV/c^2");
695 fDiff5->GetYaxis()->SetTitle("Entries");
696
697 fOutput->Add(fDiff5);
698
699 fDiffSideBand = new TH1F("DiffSide","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
700 fDiffSideBand->SetStats(kTRUE);
701 fDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
702 fDiffSideBand->GetYaxis()->SetTitle("Entries");
703
704 fOutput->Add(fDiffSideBand);
705
706 fDiffSideBand1 = new TH1F("DiffSide1","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
707 fDiffSideBand1->SetStats(kTRUE);
708 fDiffSideBand1->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
709 fDiffSideBand1->GetYaxis()->SetTitle("Entries");
710
711 fOutput->Add(fDiffSideBand1);
712
713 fDiffSideBand2 = new TH1F("DiffSide2","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
714 fDiffSideBand2->SetStats(kTRUE);
715 fDiffSideBand2->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
716 fDiffSideBand2->GetYaxis()->SetTitle("Entries");
717
718 fOutput->Add(fDiffSideBand2);
719
720 fDiffSideBand3 = new TH1F("DiffSide3","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
721 fDiffSideBand3->SetStats(kTRUE);
722 fDiffSideBand3->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
723 fDiffSideBand3->GetYaxis()->SetTitle("Entries");
724
725 fOutput->Add(fDiffSideBand3);
726
727 fDiffSideBand4 = new TH1F("DiffSide4","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
728 fDiffSideBand4->SetStats(kTRUE);
729 fDiffSideBand4->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
730 fDiffSideBand4->GetYaxis()->SetTitle("Entries");
731
732 fOutput->Add(fDiffSideBand4);
733
734 fDiffSideBand5 = new TH1F("DiffSide5","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
735 fDiffSideBand5->SetStats(kTRUE);
736 fDiffSideBand5->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
737 fDiffSideBand5->GetYaxis()->SetTitle("Entries");
738
739 fOutput->Add(fDiffSideBand5);
740
741 fDiffWrongSign = new TH1F("DiffWrong","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
742 fDiffWrongSign->SetStats(kTRUE);
743 fDiffWrongSign->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
744 fDiffWrongSign->GetYaxis()->SetTitle("Entries");
745
746 fOutput->Add(fDiffWrongSign);
747
748 fDiffWrongSign1 = new TH1F("DiffWrong1","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
749 fDiffWrongSign1->SetStats(kTRUE);
750 fDiffWrongSign1->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
751 fDiffWrongSign1->GetYaxis()->SetTitle("Entries");
752
753 fOutput->Add(fDiffWrongSign1);
754
755 fDiffWrongSign2 = new TH1F("DiffWrong2","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
756 fDiffWrongSign2->SetStats(kTRUE);
757 fDiffWrongSign2->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
758 fDiffWrongSign2->GetYaxis()->SetTitle("Entries");
759
760 fOutput->Add(fDiffWrongSign2);
761
762 fDiffWrongSign3 = new TH1F("DiffWrong3","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
763 fDiffWrongSign3->SetStats(kTRUE);
764 fDiffWrongSign3->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
765 fDiffWrongSign3->GetYaxis()->SetTitle("Entries");
766
767 fOutput->Add(fDiffWrongSign3);
768
769 fDiffWrongSign4 = new TH1F("DiffWrong4","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
770 fDiffWrongSign4->SetStats(kTRUE);
771 fDiffWrongSign4->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
772 fDiffWrongSign4->GetYaxis()->SetTitle("Entries");
773
774 fOutput->Add(fDiffWrongSign4);
775
776 fDiffWrongSign5 = new TH1F("DiffWrong5","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
777 fDiffWrongSign5->SetStats(kTRUE);
778 fDiffWrongSign5->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
779 fDiffWrongSign5->GetYaxis()->SetTitle("Entries");
780
781 fOutput->Add(fDiffWrongSign5);
782
783 fDStarMass = new TH1F("RECODStar2","RECO DStar invariant mass distribution",750,1.5,2.5);
784
785 fOutput->Add(fDStarMass);
786
787 fTrueDiff = new TH1F("dstar","True Reco diff",750,0,0.2);
788
789 fOutput->Add(fTrueDiff);
790 fTrueDiff2 = new TH2F("DiffDstar_pt","True Reco diff vs pt",100,0,15,900,0,0.3);
791
792 fOutput->Add(fTrueDiff2);
793
794 fPtDStar = new TH1F("DStarpt","Reconstructed D* candidates momentum ",500,0,25);
795 fPtDStar->SetStats(kTRUE);
796 fPtDStar->GetXaxis()->SetTitle("GeV/c");
797 fPtDStar->GetYaxis()->SetTitle("Entries");
798
799 fOutput->Add(fPtDStar);
800
801 fMCDStarPt = new TH1F("DStarMC2","Reconstructed D* momentum MC tagged ",500,0,25);
802 fMCDStarPt->SetStats(kTRUE);
803 fMCDStarPt->GetXaxis()->SetTitle("GeV/c");
804 fMCDStarPt->GetYaxis()->SetTitle("Entries");
805
806 fOutput->Add(fMCDStarPt);
807
808 return kTRUE;
809
810}
811//______________________________ side band background for D*___________________________________
812void AliAnalysisTaskSEDStarSpectra::SideBandBackground(Double_t finvM, Double_t finvMDStar, Double_t pt, Int_t okD0, Int_t okD0bar ){
813
814 // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
815 // (expected detector resolution) on the left and right frm the D0 mass. Each band
816 // has a width of ~5 sigmas. Two band needed for opening angle considerations
817
818 if((finvM>=1.740 && finvM<=1.829) || (finvM>=1.901 && finvM<=1.990)){
819
820 if(okD0==1) fDiffSideBand->Fill(finvMDStar-finvM); // M(Kpipi)-M(Kpi) side band background
821 if(okD0bar==1) fDiffSideBand->Fill(finvMDStar-finvM); // M(Kpipi)-M(Kpi) side band background
822
823 // pt bins
824 if( pt <= 1){
825 if(okD0==1) fDiffSideBand1->Fill(finvMDStar-finvM);
826 if(okD0bar==1) fDiffSideBand1->Fill(finvMDStar-finvM);
827 }
828 if( pt > 1 && pt <=3){
829 if(okD0==1 ) fDiffSideBand2->Fill(finvMDStar-finvM);
830 if(okD0bar==1 ) fDiffSideBand2->Fill(finvMDStar-finvM);
831 }
832 if( pt >3 && pt <=5){
833 if(okD0==1) fDiffSideBand3->Fill(finvMDStar-finvM);
834 if(okD0bar==1) fDiffSideBand3->Fill(finvMDStar-finvM);
835 }
836 if( pt > 5 && pt <8){
837 if(okD0==1) fDiffSideBand4->Fill(finvMDStar-finvM);
838 if(okD0bar==1) fDiffSideBand4->Fill(finvMDStar-finvM);
839 }
840 if( pt >= 8){
841 if(okD0==1) fDiffSideBand5->Fill(finvMDStar-finvM);
842 if(okD0bar==1) fDiffSideBand5->Fill(finvMDStar-finvM);
843 }
844
845 }
846}
847
848//________________________________________________________________________________________________________________
849void AliAnalysisTaskSEDStarSpectra::WrongSignForDStar(Double_t finvM, Double_t finvMDStar, Double_t pt, Int_t okD0, Int_t okD0bar){
850 //
851 // assign the wrong charge to the soft pion to create background
852 //
853 if(okD0==1 ) fDiffWrongSign->Fill(finvMDStar-finvM); // M(Kpipi)-M(Kpi) side band background
854 if(okD0bar==1 ) fDiffWrongSign->Fill(finvMDStar-finvM); // M(Kpipi)-M(Kpi) side band background
855
856 // pt bins
857 if( pt <= 1){
858 if(okD0==1) fDiffWrongSign1->Fill(finvMDStar-finvM);
859 if(okD0bar==1) fDiffWrongSign1->Fill(finvMDStar-finvM);
860 }
861 if( pt > 1 && pt <=3){
862 if(okD0==1) fDiffWrongSign2->Fill(finvMDStar-finvM);
863 if(okD0bar==1) fDiffWrongSign2->Fill(finvMDStar-finvM);
864 }
865 if( pt >3 && pt <=5){
866 if(okD0==1) fDiffWrongSign3->Fill(finvMDStar-finvM);
867 if(okD0bar==1) fDiffWrongSign3->Fill(finvMDStar-finvM);
868 }
869 if( pt > 5 && pt <8){
870 if(okD0==1) fDiffWrongSign4->Fill(finvMDStar-finvM);
871 if(okD0bar==1) fDiffWrongSign4->Fill(finvMDStar-finvM);
872 }
873 if( pt >= 8){
874 if(okD0==1) fDiffWrongSign5->Fill(finvMDStar-finvM);
875 if(okD0bar==1) fDiffWrongSign5->Fill(finvMDStar-finvM);
876 }
877
878}
879//_____________________________________________________________________________________________
fca7f2a2 880Bool_t AliAnalysisTaskSEDStarSpectra::SetUtrechtSelections(Double_t ptD0){
645e004b 881
882 //cuts[0] = inv. mass half width [GeV]
883 //cuts[1] = dca [cm]
884 //cuts[2] = cosThetaStar
885 //cuts[3] = pTK [GeV/c]
886 //cuts[4] = pTPi [GeV/c]
887 //cuts[5] = d0K [cm] upper limit!
888 //cuts[6] = d0Pi [cm] upper limit!
889 //cuts[7] = d0d0 [cm^2]
890 //cuts[8] = cosThetaPoint
891
892 if(ptD0<1.){
893 fVHF->SetD0toKpiCuts(0.450,0.04,0.8,0.21,0.21,0.021,0.021,-0.0002,0.9);
894 }
895 else if(ptD0 >=1. && ptD0 <=2.){
896 fVHF->SetD0toKpiCuts(0.450,0.02,0.7,0.8,0.8,0.021,0.021,-0.0002,0.9);
897 }
898 else if(ptD0 >2. && ptD0 <=3.){
899 fVHF->SetD0toKpiCuts(0.450,0.04,0.8,0.8,0.8,0.035,0.042,-0.000085,0.9);
900 }
901 else if(ptD0 >3. && ptD0 <=5.){
902 fVHF->SetD0toKpiCuts(0.450,0.016,0.8,1.2,1.2,0.042,0.056,-0.000085,0.9);
903 }
904 else if(ptD0 >5.){
905 fVHF->SetD0toKpiCuts(0.450,0.08,1.0,1.2,1.2,0.07,0.07,0.0001,0.9);
906 }
907 return kTRUE;
908}
909//_____________________________ pid _______________________________________-
910Bool_t AliAnalysisTaskSEDStarSpectra::SelectPID(AliAODTrack *track, Double_t nsig){//pid - K = 3
911 //TPC
912 Bool_t isKaon=kTRUE;
913 if ((track->GetStatus()&AliESDtrack::kTPCpid )==0) return isKaon;
914 AliAODPid *pid = track->GetDetPid();
915 static AliTPCPIDResponse theTPCpid;
916 Double_t nsigma = theTPCpid.GetNumberOfSigmas(track->P(),pid->GetTPCsignal(),track->GetTPCClusterMap().CountBits(), (AliPID::kKaon));
917 if (TMath::Abs(nsigma)>nsig) isKaon=kFALSE;
918 //ITS
919 //
920 //
921 //
922
923 return isKaon;
924}
925