]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/FlavourJetTasks/AliAnalysisTaskFlavourJetCorrelations.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixed a number of warnings
[u/mrichter/AliRoot.git] / PWGJE / FlavourJetTasks / AliAnalysisTaskFlavourJetCorrelations.cxx
CommitLineData
931f8b00 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// Analysis Taks for reconstructed particle correlation
18// (first implementation done for D mesons) with jets
19// (use the so called Emcal framework)
20//
21//-----------------------------------------------------------------------
22// Authors:
23// C. Bianchin (Utrecht University) chiara.bianchin@cern.ch
24// A. Grelli (Utrecht University) a.grelli@uu.nl
25// X. Zhang (LBNL) XMZhang@lbl.gov
26//-----------------------------------------------------------------------
27
28#include <TDatabasePDG.h>
29#include <TParticle.h>
30#include "TROOT.h"
31#include <TH3F.h>
32#include <THnSparse.h>
33
34#include "AliAnalysisTaskFlavourJetCorrelations.h"
35#include "AliAODMCHeader.h"
36#include "AliAODHandler.h"
37#include "AliAnalysisManager.h"
38#include "AliLog.h"
39#include "AliEmcalJet.h"
40#include "AliJetContainer.h"
41#include "AliAODRecoDecay.h"
42#include "AliAODRecoCascadeHF.h"
43#include "AliAODRecoDecayHF2Prong.h"
44#include "AliESDtrack.h"
45#include "AliAODMCParticle.h"
46#include "AliPicoTrack.h"
47#include "AliRDHFCutsD0toKpi.h"
48#include "AliRDHFCutsDStartoKpipi.h"
49
50ClassImp(AliAnalysisTaskFlavourJetCorrelations)
51
52
53//_______________________________________________________________________________
54
55AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations() :
56AliAnalysisTaskEmcalJet("",kFALSE),
57fUseMCInfo(kTRUE),
58fUseReco(kTRUE),
59fCandidateType(),
60fPDGmother(),
61fNProngs(),
62fPDGdaughters(),
63fBranchName(),
64fmyOutput(0),
65fCuts(0),
66fMinMass(),
67fMaxMass(),
68fJetArrName(0),
69fCandArrName(0),
70fLeadingJetOnly(kFALSE),
bbb94467 71fJetRadius(0),
72fCandidateArray(0),
73fSideBandArray(0)
931f8b00 74{
75 //
76 // Default ctor
77 //
78 //SetMakeGeneralHistograms(kTRUE);
79
80}
81
82//_______________________________________________________________________________
83
84AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations(const Char_t* name, AliRDHFCuts* cuts,ECandidateType candtype) :
85AliAnalysisTaskEmcalJet(name,kFALSE),
86fUseMCInfo(kTRUE),
87fUseReco(kTRUE),
88fCandidateType(),
89fPDGmother(),
90fNProngs(),
91fPDGdaughters(),
92fBranchName(),
93fmyOutput(0),
94fCuts(0),
95fMinMass(),
96fMaxMass(),
97fJetArrName(0),
98fCandArrName(0),
99fLeadingJetOnly(kFALSE),
bbb94467 100fJetRadius(0),
101fCandidateArray(0),
102fSideBandArray(0)
931f8b00 103{
104 //
105 // Constructor. Initialization of Inputs and Outputs
106 //
107
108 Info("AliAnalysisTaskFlavourJetCorrelations","Calling Constructor");
109 fCuts=cuts;
110 fCandidateType=candtype;
111 const Int_t nptbins=fCuts->GetNPtBins();
112 Float_t defaultSigmaD013[13]={0.012, 0.012, 0.012, 0.015, 0.015,0.018,0.018,0.020,0.020,0.030,0.030,0.037,0.040};
113
114 switch(fCandidateType){
115 case 0:
116 fPDGmother=421;
117 fNProngs=2;
118 fPDGdaughters[0]=211;//pi
119 fPDGdaughters[1]=321;//K
120 fPDGdaughters[2]=0; //empty
121 fPDGdaughters[3]=0; //empty
122 fBranchName="D0toKpi";
123 fCandArrName="D0";
124 break;
125 case 1:
126 fPDGmother=413;
127 fNProngs=3;
128 fPDGdaughters[1]=211;//pi soft
129 fPDGdaughters[0]=421;//D0
130 fPDGdaughters[2]=211;//pi fromD0
131 fPDGdaughters[3]=321; // K from D0
132 fBranchName="Dstar";
133 fCandArrName="Dstar";
134
135 if(nptbins<=13){
136 for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];
137 } else {
138 AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
139 }
140 break;
141 default:
142 printf("%d not accepted!!\n",fCandidateType);
143 break;
144 }
145
146 if(fCandidateType==kD0toKpi)SetMassLimits(0.15,fPDGmother);
147 if(fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
148
bbb94467 149 DefineInput(1, TClonesArray::Class());
150 DefineInput(2, TClonesArray::Class());
151
931f8b00 152 DefineOutput(1,TList::Class()); // histos
153 DefineOutput(2,AliRDHFCuts::Class()); // my cuts
154
155}
156
157//_______________________________________________________________________________
158
159AliAnalysisTaskFlavourJetCorrelations::~AliAnalysisTaskFlavourJetCorrelations() {
160 //
161 // destructor
162 //
163
164 Info("~AliAnalysisTaskFlavourJetCorrelations","Calling Destructor");
165
166 delete fmyOutput;
167 delete fCuts;
168
169}
170
171//_______________________________________________________________________________
172
173void AliAnalysisTaskFlavourJetCorrelations::Init(){
174 //
175 // Initialization
176 //
177
178 if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");
179 switch(fCandidateType){
180 case 0:
181 {
182 AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
183 copyfCuts->SetName("AnalysisCutsDzero");
184 // Post the data
185 PostData(2,copyfCuts);
186 }
187 break;
188 case 1:
189 {
190 AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
191 copyfCuts->SetName("AnalysisCutsDStar");
192 // Post the cuts
193 PostData(2,copyfCuts);
194 }
195 break;
196 default:
197 return;
198 }
199
200 return;
201}
202
203//_______________________________________________________________________________
204
205void AliAnalysisTaskFlavourJetCorrelations::UserCreateOutputObjects() {
206 // output
207 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
208 fmyOutput = new TList();
209 fmyOutput->SetOwner();
210 fmyOutput->SetName("pippo");
211 // define histograms
212 DefineHistoForAnalysis();
213 PostData(1,fmyOutput);
214
215 return;
216}
217
218//_______________________________________________________________________________
219
220void AliAnalysisTaskFlavourJetCorrelations::UserExec(Option_t *)
221{
222 // user exec
223 if (!fInitialized){
224 AliAnalysisTaskEmcalJet::ExecOnce();
225 }
226
227 // Load the event
228 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
229
230 TClonesArray *arrayDStartoD0pi=0;
231 TClonesArray *trackArr = 0;
931f8b00 232
233 if (!aodEvent && AODEvent() && IsStandardAOD()) {
234
235 // In case there is an AOD handler writing a standard AOD, use the AOD
236 // event in memory rather than the input (ESD) event.
237 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
238
239 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
240 // have to taken from the AOD event hold by the AliAODExtension
241 AliAODHandler* aodHandler = (AliAODHandler*)
242 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
243 if(aodHandler->GetExtensions()) {
244 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
245 AliAODEvent *aodFromExt = ext->GetAOD();
246 arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
247 }
248 } else {
249 arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
250 }
251
252 if (!arrayDStartoD0pi) {
253 AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
254 // return;
255 } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
256
257 TClonesArray* mcArray = 0x0;
258 if (fUseMCInfo) {
259 mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
260 if (!mcArray) {
261 printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
262 return;
263 }
264 }
265
266 //retrieve jets
267 trackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoTracks"));
268 //clusArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClustersCorr"));
269 //jetArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName));
270 fJetRadius=GetJetContainer(0)->GetJetRadius();
271
272 if(!trackArr){
273 AliInfo(Form("Could not find the track array\n"));
274 return;
275 }
276
bbb94467 277
c3c1ab31 278 fCandidateArray = dynamic_cast<TClonesArray*>(GetInputData(1));
279 if (!fCandidateArray) return;
280 if (fCandidateType==1) {
281 fSideBandArray = dynamic_cast<TClonesArray*>(GetInputData(2));
282 if (!fSideBandArray) return;
283 }
284 //Printf("ncandidates found %d",fCandidateArray->GetEntriesFast());
931f8b00 285
286 //Histograms
287 TH1I* hstat=(TH1I*)fmyOutput->FindObject("hstat");
288 TH1F* hPtJetTrks=(TH1F*)fmyOutput->FindObject("hPtJetTrks");
289 TH1F* hPhiJetTrks=(TH1F*)fmyOutput->FindObject("hPhiJetTrks");
290 TH1F* hEtaJetTrks=(TH1F*)fmyOutput->FindObject("hEtaJetTrks");
291 TH1F* hEjetTrks=(TH1F*)fmyOutput->FindObject("hEjetTrks");
292 TH1F* hPtJet=(TH1F*)fmyOutput->FindObject("hPtJet");
293 TH1F* hPhiJet=(TH1F*)fmyOutput->FindObject("hPhiJet");
294 TH1F* hEtaJet=(TH1F*)fmyOutput->FindObject("hEtaJet");
295 TH1F* hEjet=(TH1F*)fmyOutput->FindObject("hEjet");
296 TH1F* hNtrArr=(TH1F*)fmyOutput->FindObject("hNtrArr");
297 TH1F* hNJetPerEv=(TH1F*)fmyOutput->FindObject("hNJetPerEv");
d9126a4c
CB		 298 TH1F* hdeltaRJetTracks=(TH1F*)fmyOutput->FindObject("hdeltaRJetTracks");
299 TH1F* hNDPerEvNoJet=(TH1F*)fmyOutput->FindObject("hNDPerEvNoJet");
300 TH1F* hptDPerEvNoJet=(TH1F*)fmyOutput->FindObject("hptDPerEvNoJet");
301 TH1F* hNJetPerEvNoD=(TH1F*)fmyOutput->FindObject("hNJetPerEvNoD");
1353fa73 302 TH1F* hPtJetPerEvNoD=(TH1F*)fmyOutput->FindObject("hPtJetPerEvNoD");
303 THnSparseF* hnspDstandalone=(THnSparseF*)fmyOutput->FindObject("hsDstandalone");
d9126a4c 304
931f8b00 305 hstat->Fill(0);
306
307 // fix for temporary bug in ESDfilter
308 // the AODs with null vertex pointer didn't pass the PhysSel
309 if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
310
311 //Event selection
312 Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
313 TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
314 if(!iseventselected) return;
315
316 hstat->Fill(1);
931f8b00 317
318 //retrieve charm candidates selected
bbb94467 319 Int_t candidates = fCandidateArray->GetEntriesFast();
d9126a4c
CB		 320
321 //trigger on jets
931f8b00 322
d9126a4c
CB		 323 Int_t njets=GetJetContainer()->GetNJets();
324 if(njets == 0) {
325 hstat->Fill(6, candidates);
326 hNDPerEvNoJet->Fill(candidates);
327 for(Int_t iD=0;iD<candidates;iD++){
bbb94467 328 AliVParticle* cand=(AliVParticle*)fCandidateArray->At(iD);
d9126a4c
CB		 329 if(!cand) continue;
330 hptDPerEvNoJet->Fill(cand->Pt());
331
332 }
333 return;
334
335 }
1353fa73 336
337 //loop on candidates standalone (checking the candidates are there and their phi-eta distributions)
338
339 for(Int_t ic = 0; ic < candidates; ic++) {
340
341 // D* candidates
bbb94467 342 AliAODRecoDecayHF* charm=0x0;
343 AliAODRecoCascadeHF* dstar=0x0;
344
345
346 charm=(AliAODRecoDecayHF*)fCandidateArray->At(ic);
1353fa73 347 if(!charm) continue;
bbb94467 348
349 if(fCandidateType==kDstartoKpipi){
350 dstar=(AliAODRecoCascadeHF*)fCandidateArray->At(ic);
351 if(!dstar) continue;
352 }
353
1353fa73 354 hstat->Fill(2);
355
356 Double_t candsparse[4]={charm->Eta(), charm->Phi(), charm->Pt(), 0};
357
358 if(fCandidateType==kDstartoKpipi) {
1353fa73 359 Double_t deltamass= dstar->DeltaInvMass();
360 candsparse[3]=deltamass;
361 hnspDstandalone->Fill(candsparse);
362 }
363 if(fCandidateType==kD0toKpi){
364 AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)charm;
365 Int_t isselected=fCuts->IsSelected(dzero,AliRDHFCuts::kAll,aodEvent);
366
367 Double_t masses[2];
368 Int_t pdgdaughtersD0[2]={211,321};//pi,K
369 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
370
371 masses[0]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
372 masses[1]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
373 if(isselected==1 || isselected==3) {
374 candsparse[3]=masses[0];
375 hnspDstandalone->Fill(candsparse);
376 }
377 if(isselected>=2){
378 candsparse[3]=masses[1];
379 hnspDstandalone->Fill(candsparse);
380
381 }
382 }
383 }
384
931f8b00 385 // we start with jets
386 Double_t ejet = 0;
387 Double_t phiJet = 0;
388 Double_t etaJet = 0;
389 Double_t ptjet = 0;
390 Double_t leadingJet =0;
391
392 Int_t ntrarr=trackArr->GetEntriesFast();
393 hNtrArr->Fill(ntrarr);
394
395 for(Int_t i=0;i<ntrarr;i++){
396 AliVTrack *jtrack=static_cast<AliVTrack*>(trackArr->At(i));
397 //reusing histograms
398 hPtJetTrks->Fill(jtrack->Pt());
399 hPhiJetTrks->Fill(jtrack->Phi());
400 hEtaJetTrks->Fill(jtrack->Eta());
401 hEjetTrks->Fill(jtrack->E());
402 }
403
404
405 //option to use only the leading jet
406 if(fLeadingJetOnly){
407 for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {
408 AliEmcalJet* jetL = (AliEmcalJet*)GetJetFromArray(iJetsL);
409 ptjet = jetL->Pt();
410 if(ptjet>leadingJet ) leadingJet = ptjet;
411
412 }
413 }
414
415 Int_t cntjet=0;
416 //loop over jets and consider only the leading jet in the event
417 for (Int_t iJets = 0; iJets<njets; iJets++) {
418 //Printf("Jet N %d",iJets);
419 AliEmcalJet* jet = (AliEmcalJet*)GetJetFromArray(iJets);
420 //Printf("Corr task Accept Jet");
421 if(!AcceptJet(jet)) {
422 hstat->Fill(5);
423 continue;
424 }
425
426 //jets variables
427 ejet = jet->E();
428 phiJet = jet->Phi();
429 etaJet = jet->Eta();
430 ptjet = jet->Pt();
431
432 // choose the leading jet
433 if(fLeadingJetOnly && (ejet<leadingJet)) continue;
434 //used for normalization
435 hstat->Fill(3);
436 cntjet++;
437 // fill energy, eta and phi of the jet
438 hEjet ->Fill(ejet);
439 hPhiJet ->Fill(phiJet);
440 hEtaJet ->Fill(etaJet);
441 hPtJet ->Fill(ptjet);
442
443 //loop on jet particles
444 Int_t ntrjet= jet->GetNumberOfTracks();
445 for(Int_t itrk=0;itrk<ntrjet;itrk++){
446
447 AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,trackArr);
448 hdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));
449
450 }//end loop on jet tracks
451
d9126a4c
CB		 452 if(candidates==0){
453 hstat->Fill(7);
454 hPtJetPerEvNoD->Fill(jet->Pt());
455 }
456
931f8b00 457 //Printf("N candidates %d ", candidates);
458 for(Int_t ic = 0; ic < candidates; ic++) {
459
460 // D* candidates
461 AliVParticle* charm=0x0;
bbb94467 462 charm=(AliVParticle*)fCandidateArray->At(ic);
931f8b00 463 if(!charm) continue;
1353fa73 464
931f8b00 465
466 FlagFlavour(charm, jet);
467 if (jet->TestFlavourTag(AliEmcalJet::kDStar)) hstat->Fill(4);
468
469 FillHistogramsRecoJetCorr(charm, jet);
470
471 }
472 //retrieve side band background candidates for Dstar
bbb94467 473 Int_t nsbcand = 0; if(fSideBandArray) nsbcand=fSideBandArray->GetEntriesFast();
931f8b00 474
475 for(Int_t ib=0;ib<nsbcand;ib++){
476 if(fCandidateType==kDstartoKpipi && !fUseMCInfo){
bbb94467 477 AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)fSideBandArray->At(ib);
931f8b00 478 if(!sbcand) continue;
479 SideBandBackground(sbcand,jet);
480 }
481 if(fUseMCInfo){
482 AliAODRecoDecayHF* charmbg = 0x0;
bbb94467 483 charmbg=(AliAODRecoDecayHF*)fCandidateArray->At(ib);
931f8b00 484 if(!charmbg) continue;
485 MCBackground(charmbg,jet);
486 }
487 }
488 } // end of jet loop
489
490 hNJetPerEv->Fill(cntjet);
d9126a4c 491 if(candidates==0) hNJetPerEvNoD->Fill(cntjet);
931f8b00 492 PostData(1,fmyOutput);
493
494}
495
496//_______________________________________________________________________________
497
498void AliAnalysisTaskFlavourJetCorrelations::Terminate(Option_t*)
499{
500 // The Terminate() function is the last function to be called during
501 // a query. It always runs on the client, it can be used to present
502 // the results graphically or save the results to file.
503
504 Info("Terminate"," terminate");
505 AliAnalysisTaskSE::Terminate();
506
507 fmyOutput = dynamic_cast<TList*> (GetOutputData(1));
508 if (!fmyOutput) {
509 printf("ERROR: fmyOutput not available\n");
510 return;
511 }
512}
513
514//_______________________________________________________________________________
515
516void AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){
517 Float_t mass=0;
518 Int_t abspdg=TMath::Abs(pdg);
519
520 mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
521 // compute the Delta mass for the D*
522 if(fCandidateType==kDstartoKpipi){
523 Float_t mass1=0;
524 mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
525 mass = mass-mass1;
526 }
527
528 fMinMass = mass-range;
529 fMaxMass = mass+range;
530
531 AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
532 if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
533}
534
535//_______________________________________________________________________________
536
537void AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){
538 if(uplimit>lowlimit)
539 {
540 fMinMass = lowlimit;
541 fMaxMass = uplimit;
542 }
543 else{
544 printf("Error! Lower limit larger than upper limit!\n");
545 fMinMass = uplimit - uplimit*0.2;
546 fMaxMass = uplimit;
547 }
548}
549
550//_______________________________________________________________________________
551
552Bool_t AliAnalysisTaskFlavourJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){
553 if(nptbins>30) {
554 AliInfo("Maximum number of bins allowed is 30!");
555 return kFALSE;
556 }
557 if(!width) return kFALSE;
558 for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
559 return kTRUE;
560}
561
562//_______________________________________________________________________________
563
564Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet) const{
565 if(!part || !jet){
566 printf("Particle or jet do not exist!\n");
567 return -99;
568 }
569 Double_t p[3],pj[3];
570 Bool_t okpp=part->PxPyPz(p);
571 Bool_t okpj=jet->PxPyPz(pj);
572 if(!okpp || !okpj){
573 printf("Problems getting momenta\n");
574 return -999;
575 }
576 Double_t pjet=jet->P();
577 Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet*pjet);
578 return z;
579}
580
581//_______________________________________________________________________________
582
583Bool_t AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
584
585 // Statistics
d9126a4c 586 TH1I* hstat=new TH1I("hstat","Statistics",8,-0.5,7.5);
931f8b00 587 hstat->GetXaxis()->SetBinLabel(1,"N ev anal");
588 hstat->GetXaxis()->SetBinLabel(2,"N ev sel");
d9126a4c 589 hstat->GetXaxis()->SetBinLabel(3,"N cand sel & jet");
931f8b00 590 hstat->GetXaxis()->SetBinLabel(4,"N jets");
591 hstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
592 hstat->GetXaxis()->SetBinLabel(6,"N jet rej");
d9126a4c
CB		 593 hstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
594 hstat->GetXaxis()->SetBinLabel(8,"N jets & !D");
931f8b00 595 hstat->SetNdivisions(1);
596 fmyOutput->Add(hstat);
597
598 const Int_t nbinsmass=200;
599 const Int_t nbinsptjet=500;
600 const Int_t nbinsptD=100;
601 const Int_t nbinsz=100;
602 const Int_t nbinsphi=200;
1353fa73 603 const Int_t nbinseta=100;
931f8b00 604
605 const Float_t ptjetlims[2]={0.,200.};
606 const Float_t ptDlims[2]={0.,50.};
607 const Float_t zlims[2]={0.,1.2};
608 const Float_t philims[2]={0.,6.3};
1353fa73 609 const Float_t etalims[2]={-1.5,1.5};
931f8b00 610
611 if(fCandidateType==kDstartoKpipi)
612 {
613
614 TH2F* hDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
615 hDiffSideBand->SetStats(kTRUE);
616 hDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
617 hDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
618 hDiffSideBand->Sumw2();
619 fmyOutput->Add(hDiffSideBand);
620
621
622 TH1F* hPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
623 hPtPion->SetStats(kTRUE);
624 hPtPion->GetXaxis()->SetTitle("GeV/c");
625 hPtPion->GetYaxis()->SetTitle("Entries");
626 hPtPion->Sumw2();
627 fmyOutput->Add(hPtPion);
628
629 }
d9126a4c 630 // D related histograms
1353fa73 631 TH1F *hNDPerEvNoJet=new TH1F("hNDPerEvNoJet","Number of candidates per event with no jets; N candidate/ev with no jet", 20, 0., 20.);
632 hNDPerEvNoJet->Sumw2();
633 fmyOutput->Add(hNDPerEvNoJet);
634
635 TH1F *hptDPerEvNoJet=new TH1F("hptDPerEvNoJet","pt distribution of candidates per events with no jets; p_{t}^{D} (GeV/c)",nbinsptD, ptDlims[0],ptDlims[1]);
636 hptDPerEvNoJet->Sumw2();
637 fmyOutput->Add(hptDPerEvNoJet);
638
639 const Int_t nAxisD=4;
e3017b3f 640 const Int_t nbinsSparseD[nAxisD]={nbinseta,nbinsphi,nbinsptD,nbinsmass};
641 const Double_t minSparseD[nAxisD] ={etalims[0],philims[0],ptDlims[0],fMinMass};
642 const Double_t maxSparseD[nAxisD] ={etalims[1],philims[1],ptDlims[1],fMaxMass};
1353fa73 643 THnSparseF *hsDstandalone=new THnSparseF("hsDstandalone","#phi, #eta, p_{T}^{D}, and mass", nAxisD, nbinsSparseD, minSparseD, maxSparseD);
644 hsDstandalone->Sumw2();
645
646 fmyOutput->Add(hsDstandalone);
647
931f8b00 648 // jet related fistograms
649
650 TH1F* hEjetTrks = new TH1F("hEjetTrks", "Jet tracks energy distribution;Energy (GeV)",500,0,200);
651 hEjetTrks->Sumw2();
652 TH1F* hPhiJetTrks = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
653 hPhiJetTrks->Sumw2();
1353fa73 654 TH1F* hEtaJetTrks = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
931f8b00 655 hEtaJetTrks->Sumw2();
656 TH1F* hPtJetTrks = new TH1F("hPtJetTrks", "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
657 hPtJetTrks->Sumw2();
658
659 TH1F* hEjet = new TH1F("hEjet", "Jet energy distribution;Energy (GeV)",500,0,200);
660 hEjet->Sumw2();
661 TH1F* hPhiJet = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
662 hPhiJet->Sumw2();
1353fa73 663 TH1F* hEtaJet = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
931f8b00 664 hEtaJet->Sumw2();
665 TH1F* hPtJet = new TH1F("hPtJet", "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
666 hPtJet->Sumw2();
667
668 TH3F* hPtJetWithD=new TH3F("hPtJetWithD","D-Jet Pt distribution; p_{T} (GeV/c);delta mass (GeV/c^{2})",nbinsptjet,ptjetlims[0],ptjetlims[1],nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
669 hPtJetWithD->Sumw2();
670 //for the MC this histogram is filled with the real background
671 TH3F* hPtJetWithDsb=new TH3F("hPtJetWithDsb","D(background)-Jet Pt distribution; p_{T} (GeV/c);delta mass (GeV/c^{2});p_{T}^{D} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1],nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
672 hPtJetWithDsb->Sumw2();
673
674 TH1F* hdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);
675 hdeltaRJetTracks->Sumw2();
676
677 TH1F* hNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);
678 hNtrArr->Sumw2();
d9126a4c 679
931f8b00 680 TH1F *hNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
681 hNJetPerEv->Sumw2();
682
d9126a4c
CB		 683 TH1F *hNJetPerEvNoD=new TH1F("hNJetPerEvNoD","Number of jets per event with no D; number of jets/ev with no D",10,-0.5,9.5);
684 hNJetPerEvNoD->Sumw2();
685
686 TH1F *hPtJetPerEvNoD=new TH1F("hPtJetPerEvNoD","pt distribution of jets per event with no D; p_{T}^{jet} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
687 hPtJetPerEvNoD->Sumw2();
688
931f8b00 689 fmyOutput->Add(hEjetTrks);
690 fmyOutput->Add(hPhiJetTrks);
691 fmyOutput->Add(hEtaJetTrks);
692 fmyOutput->Add(hPtJetTrks);
693 fmyOutput->Add(hEjet);
694 fmyOutput->Add(hPhiJet);
695 fmyOutput->Add(hEtaJet);
696 fmyOutput->Add(hPtJet);
697 fmyOutput->Add(hPtJetWithD);
698 fmyOutput->Add(hPtJetWithDsb);
699 fmyOutput->Add(hdeltaRJetTracks);
700 fmyOutput->Add(hNtrArr);
701 fmyOutput->Add(hNJetPerEv);
d9126a4c
CB		 702 fmyOutput->Add(hNJetPerEvNoD);
703 fmyOutput->Add(hPtJetPerEvNoD);
931f8b00 704
705 TH1F* hDeltaRD=new TH1F("hDeltaRD","#Delta R distribution of D candidates selected;#Delta R",200, 0.,10.);
706 hDeltaRD->Sumw2();
707 fmyOutput->Add(hDeltaRD);
708
709 //background (side bands for the Dstar and like sign for D0)
710 fJetRadius=GetJetContainer(0)->GetJetRadius();
711 TH2F* hInvMassptD = new TH2F("hInvMassptD",Form("D (Delta R < %.1f) invariant mass distribution p_{T}^{j} > threshold",fJetRadius),nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
712 hInvMassptD->SetStats(kTRUE);
713 hInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
714 hInvMassptD->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
715 hInvMassptD->Sumw2();
716
717 fmyOutput->Add(hInvMassptD);
718
719 if(fUseMCInfo){
720 TH2F* hInvMassptDbg = new TH2F("hInvMassptDbg",Form("Background D (Delta R < %.1f) invariant mass distribution p_{T}^{j} > threshold",fJetRadius),nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
721 hInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
722 hInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
723 hInvMassptDbg->Sumw2();
724 fmyOutput->Add(hInvMassptDbg);
725
726 }
727 Double_t pi=TMath::Pi(), philims2[2];
728 philims2[0]=-pi*1./2.;
729 philims2[1]=pi*3./2.;
ed5a88ff 730 const Int_t nAxis=7;
731 const Int_t nbinsSparse[nAxis]={nbinsz,nbinsphi,nbinsptjet,nbinsptD,nbinsmass,2, 2};
732 const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5, -0.5};
733 const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5};
734 THnSparseF *hsDphiz=new THnSparseF("hsDphiz","Z and #Delta#phi vs p_{T}^{jet}, p_{T}^{D}, mass. For SB or not and within the jet cone or not", nAxis, nbinsSparse, minSparse, maxSparse);
931f8b00 735 hsDphiz->Sumw2();
736
737 fmyOutput->Add(hsDphiz);
738
739 PostData(1, fmyOutput);
740
741 return kTRUE;
742}
743
744//_______________________________________________________________________________
745
746void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet){
747
748 Double_t ptD=candidate->Pt();
749 Double_t ptjet=jet->Pt();
750 Double_t deltaR=DeltaR(candidate,jet);
751 Double_t deltaphi = jet->Phi()-candidate->Phi();
752 if(deltaphi<=-(TMath::Pi())/2) deltaphi = deltaphi+2*(TMath::Pi());
753 if(deltaphi>(3*(TMath::Pi()))/2) deltaphi = deltaphi-2*(TMath::Pi());
754 Double_t z=Z(candidate,jet);
755
756 TH1F* hDeltaRD=(TH1F*)fmyOutput->FindObject("hDeltaRD");
757 hDeltaRD->Fill(deltaR);
758 if(fUseReco){
759 if(fCandidateType==kD0toKpi) {
760 AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;
761 FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,ptjet,deltaR, AODEvent());
762
763 }
764
765 if(fCandidateType==kDstartoKpipi) {
766 AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
767 FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,ptjet,deltaR);
768
769 }
770 } else{ //generated level
771 //AliInfo("Non reco");
772 FillHistogramsMCGenDJetCorr(deltaphi,z,ptD,ptjet,deltaR);
773 }
774}
775
776//_______________________________________________________________________________
777
778void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj,Double_t deltaR, AliAODEvent* aodEvent){
779
931f8b00 780
781 Double_t masses[2]={0.,0.};
782 Int_t pdgdaughtersD0[2]={211,321};//pi,K
783 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
784
785 masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
786 masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
787
788 TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");
789 THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
ed5a88ff 790 Double_t point[7]={z,dPhi,ptj,ptD,masses[0],0, deltaR<fJetRadius ? 1 : 0};
931f8b00 791
792 Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
793 if(isselected==1 || isselected==3) {
794
795 if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[0],ptD);
796
797 FillMassHistograms(masses[0], ptD, deltaR);
798 hsDphiz->Fill(point,1.);
799 }
800 if(isselected>=2) {
801 if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[1],ptD);
802
803 FillMassHistograms(masses[1], ptD, deltaR);
804 point[4]=masses[1];
805 hsDphiz->Fill(point,1.);
806 }
807
808}
809
810//_______________________________________________________________________________
811
812void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Double_t deltaR){
813 //dPhi and z not used at the moment,but will be (re)added
814
815 AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
816 Double_t deltamass= dstar->DeltaInvMass();
817 //Double_t massD0= dstar->InvMassD0();
818
819
820 TH1F* hPtPion=(TH1F*)fmyOutput->FindObject("hPtPion");
821 hPtPion->Fill(softpi->Pt());
822
823 TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");
824 THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
ed5a88ff 825 Int_t isSB=0;//IsDzeroSideBand(dstar);
931f8b00 826
ed5a88ff 827 Double_t point[7]={z,dPhi,ptj,ptD,deltamass,isSB, deltaR<fJetRadius ? 1 : 0};
931f8b00 828
829 if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,deltamass,ptD);
830
831 FillMassHistograms(deltamass, ptD, deltaR);
832 hsDphiz->Fill(point,1.);
833
834}
835
836//_______________________________________________________________________________
837
838void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t z,Double_t ptD,Double_t ptjet,Double_t deltaR){
839
840 Double_t pdgmass=0;
841 TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");
842 THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
ed5a88ff 843 Double_t point[7]={z,dPhi,ptjet,ptD,pdgmass,0, deltaR<fJetRadius ? 1 : 0};
931f8b00 844
845 if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
846 if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
847 point[4]=pdgmass;
848
849 if(deltaR<fJetRadius) {
850 hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
851 hsDphiz->Fill(point,1.);
852 }
853
854}
855
856//_______________________________________________________________________________
857
858void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD, Double_t deltaR){
859
860 if(deltaR<fJetRadius) {
861 TH2F* hInvMassptD=(TH2F*)fmyOutput->FindObject("hInvMassptD");
862 hInvMassptD->Fill(mass,ptD);
863 }
864}
865
1353fa73 866//________________________________________________________________________________
867
931f8b00 868void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliVParticle *charm, AliEmcalJet *jet){
869 Double_t deltaR=DeltaR(charm, jet);
870 AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
871 if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;
872 if (deltaR<fJetRadius) jet->AddFlavourTag(tag);
873
874 return;
875
876}
877
878//_______________________________________________________________________________
879
880void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){
881
882 // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
883 // (expected detector resolution) on the left and right frm the D0 mass. Each band
884 // has a width of ~5 sigmas. Two band needed for opening angle considerations
885 TH2F* hDiffSideBand=(TH2F*)fmyOutput->FindObject("hDiffSideBand");
886 TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb");
ed5a88ff 887 THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
931f8b00 888
889 Double_t deltaM=candDstar->DeltaInvMass();
890 //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
ed5a88ff 891 Double_t z=Z(candDstar,jet);
931f8b00 892 Double_t ptD=candDstar->Pt();
893 Double_t ptjet=jet->Pt();
894 Double_t dPhi=jet->Phi()-candDstar->Phi();
895 Double_t deltaR=DeltaR(candDstar,jet);
896 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
897 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
898
ed5a88ff 899 Int_t isSideBand=1;//IsDzeroSideBand(candDstar);
900 Double_t point[7]={z,dPhi,ptjet,ptD,deltaM,isSideBand, deltaR<fJetRadius ? 1 : 0};
931f8b00 901 //fill the background histogram with the side bands or when looking at MC with the real background
902 if(isSideBand==1){
903 hDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background
904 //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
ed5a88ff 905 hsDphiz->Fill(point,1.);
931f8b00 906 if(deltaR<fJetRadius){ // evaluate in the near side
907 //hzptDB->Fill(Z(candDstar,jet),deltaM,ptD);
908 hPtJetWithDsb->Fill(ptjet,deltaM,ptD);
909 }
910 }
911}
912
913//_______________________________________________________________________________
914
915void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg, AliEmcalJet *jet){
916
917 //need mass, deltaR, pt jet, ptD
918 //two cases: D0 and Dstar
919
920 Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());
921 if(!isselected) return;
922
923 Double_t ptD=candbg->Pt();
924 Double_t ptjet=jet->Pt();
925 Double_t deltaR=DeltaR(candbg,jet);
926
927 TH2F* hInvMassptDbg=(TH2F*)fmyOutput->FindObject("hInvMassptDbg");
928 TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb");
929
930
931 if(fCandidateType==kDstartoKpipi){
932 AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
933 Double_t deltaM=dstarbg->DeltaInvMass();
934 hInvMassptDbg->Fill(deltaM,ptD);
935 if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,deltaM,ptD);
936 }
937
938 if(fCandidateType==kD0toKpi){
939 Double_t masses[2]={0.,0.};
940 Int_t pdgdaughtersD0[2]={211,321};//pi,K
941 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
942
943 masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
944 masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
945
946
947 if(isselected==1 || isselected==3) {
948 if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[0],ptD);
949 hInvMassptDbg->Fill(masses[0],ptD);
950 }
951 if(isselected>=2) {
952 if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[1],ptD);
953 hInvMassptDbg->Fill(masses[1],ptD);
954 }
955
956
957 }
958}
959
960//_______________________________________________________________________________
961
962Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliVParticle *p1, AliVParticle *p2) const {
963 //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
964
965 if(!p1 || !p2) return -1;
966 Double_t phi1=p1->Phi(),eta1=p1->Eta();
967 Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
968
969 Double_t dPhi=phi1-phi2;
970 if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
971 if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
972
973 Double_t dEta=eta1-eta2;
974 Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
975 return deltaR;
976
977}
978
979//_______________________________________________________________________________
980
981Int_t AliAnalysisTaskFlavourJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){
982
983 Double_t ptD=candDstar->Pt();
984 Int_t bin = fCuts->PtBin(ptD);
985 if (bin < 0){
986 // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
987 bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
988 return -1; // out of bounds
989 }
990
991 Double_t invM=candDstar->InvMassD0();
992 Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
993
994 Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
995 Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
996
997 if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
998 else return 0;
999
1000}