]>
Commit | Line | Data |
---|---|---|
ae040693 | 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 "AliAnalysisTaskFlavourFilterAndJetCorrelations.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 | ||
50 | ClassImp(AliAnalysisTaskFlavourFilterAndJetCorrelations) | |
51 | ||
52 | ||
53 | //_______________________________________________________________________________ | |
54 | ||
55 | AliAnalysisTaskFlavourFilterAndJetCorrelations::AliAnalysisTaskFlavourFilterAndJetCorrelations() : | |
56 | AliAnalysisTaskEmcalJet("",kFALSE), | |
57 | fUseMCInfo(kTRUE), | |
58 | fUseReco(kTRUE), | |
59 | fCandidateType(), | |
60 | fPDGmother(), | |
61 | fNProngs(), | |
62 | fPDGdaughters(), | |
63 | fBranchName(), | |
64 | fmyOutput(0), | |
963e733a | 65 | fmyOutputF(0), |
ae040693 | 66 | fCuts(0), |
67 | fMinMass(), | |
68 | fMaxMass(), | |
69 | fJetArrName(0), | |
70 | fCandArrName(0), | |
71 | fLeadingJetOnly(kFALSE), | |
963e733a | 72 | fJetRadius(0), |
73 | fCandidateArray(0), | |
74 | fSideBandArray(0) | |
ae040693 | 75 | { |
76 | // | |
77 | // Default ctor | |
78 | // | |
79 | //SetMakeGeneralHistograms(kTRUE); | |
80 | ||
81 | } | |
82 | ||
83 | //_______________________________________________________________________________ | |
84 | ||
85 | AliAnalysisTaskFlavourFilterAndJetCorrelations::AliAnalysisTaskFlavourFilterAndJetCorrelations(const Char_t* name, AliRDHFCuts* cuts,ECandidateType candtype) : | |
86 | AliAnalysisTaskEmcalJet(name,kFALSE), | |
87 | fUseMCInfo(kTRUE), | |
88 | fUseReco(kTRUE), | |
89 | fCandidateType(), | |
90 | fPDGmother(), | |
91 | fNProngs(), | |
92 | fPDGdaughters(), | |
93 | fBranchName(), | |
94 | fmyOutput(0), | |
95 | fmyOutputF(0), | |
96 | fCuts(0), | |
97 | fMinMass(), | |
98 | fMaxMass(), | |
99 | fJetArrName(0), | |
100 | fCandArrName(0), | |
101 | fLeadingJetOnly(kFALSE), | |
963e733a | 102 | fJetRadius(0), |
103 | fCandidateArray(0), | |
104 | fSideBandArray(0) | |
ae040693 | 105 | { |
106 | // | |
107 | // Constructor. Initialization of Inputs and Outputs | |
108 | // | |
109 | ||
110 | Info("AliAnalysisTaskFlavourFilterAndJetCorrelations","Calling Constructor"); | |
111 | fCuts=cuts; | |
112 | fCandidateType=candtype; | |
113 | const Int_t nptbins=fCuts->GetNPtBins(); | |
114 | 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}; | |
115 | ||
116 | switch(fCandidateType){ | |
117 | case 0: | |
118 | fPDGmother=421; | |
119 | fNProngs=2; | |
120 | fPDGdaughters[0]=211;//pi | |
121 | fPDGdaughters[1]=321;//K | |
122 | fPDGdaughters[2]=0; //empty | |
123 | fPDGdaughters[3]=0; //empty | |
124 | fBranchName="D0toKpi"; | |
125 | fCandArrName="D0"; | |
126 | break; | |
127 | case 1: | |
128 | fPDGmother=413; | |
129 | fNProngs=3; | |
130 | fPDGdaughters[1]=211;//pi soft | |
131 | fPDGdaughters[0]=421;//D0 | |
132 | fPDGdaughters[2]=211;//pi fromD0 | |
133 | fPDGdaughters[3]=321; // K from D0 | |
134 | fBranchName="Dstar"; | |
135 | fCandArrName="Dstar"; | |
136 | ||
137 | if(nptbins<=13){ | |
138 | for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt]; | |
139 | } else { | |
140 | AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins)); | |
141 | } | |
142 | break; | |
143 | default: | |
144 | printf("%d not accepted!!\n",fCandidateType); | |
145 | break; | |
146 | } | |
147 | ||
148 | if(fCandidateType==kD0toKpi)SetMassLimits(0.15,fPDGmother); | |
149 | if(fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother); | |
150 | ||
151 | DefineOutput(1,TList::Class()); // histos | |
152 | DefineOutput(2,AliRDHFCuts::Class()); // my cuts | |
153 | DefineOutput(3,TList::Class()); //histos filter | |
154 | ||
155 | } | |
156 | ||
157 | //_______________________________________________________________________________ | |
158 | ||
159 | AliAnalysisTaskFlavourFilterAndJetCorrelations::~AliAnalysisTaskFlavourFilterAndJetCorrelations() { | |
160 | // | |
161 | // destructor | |
162 | // | |
163 | ||
164 | Info("~AliAnalysisTaskFlavourFilterAndJetCorrelations","Calling Destructor"); | |
165 | ||
166 | delete fmyOutput; | |
167 | delete fmyOutputF; | |
168 | delete fCuts; | |
169 | ||
170 | } | |
171 | ||
172 | //_______________________________________________________________________________ | |
173 | ||
174 | void AliAnalysisTaskFlavourFilterAndJetCorrelations::Init(){ | |
175 | // | |
176 | // Initialization | |
177 | // | |
178 | ||
179 | if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n"); | |
180 | switch(fCandidateType){ | |
181 | case 0: | |
182 | { | |
183 | AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts))); | |
184 | copyfCuts->SetName("AnalysisCutsDzero"); | |
185 | // Post the data | |
186 | PostData(2,copyfCuts); | |
187 | } | |
188 | break; | |
189 | case 1: | |
190 | { | |
191 | AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts))); | |
192 | copyfCuts->SetName("AnalysisCutsDStar"); | |
193 | // Post the cuts | |
194 | PostData(2,copyfCuts); | |
195 | } | |
196 | break; | |
197 | default: | |
198 | return; | |
199 | } | |
200 | ||
201 | return; | |
202 | } | |
203 | ||
204 | //_______________________________________________________________________________ | |
205 | ||
206 | void AliAnalysisTaskFlavourFilterAndJetCorrelations::UserCreateOutputObjects() { | |
963e733a | 207 | |
208 | fCandidateArray = new TClonesArray("AliAODRecoDecayHF",0); | |
209 | fCandidateArray->SetName(Form("fCandidateArray%s%s",fCandArrName.Data(),fUseReco ? "rec" : "gen")); | |
210 | ||
211 | if (fCandidateType==kDstartoKpipi){ | |
212 | fSideBandArray = new TClonesArray("AliAODRecoCascadeHF",0); //this is for the DStar only! | |
213 | fSideBandArray->SetName(Form("fSideBandArray%s%s",fCandArrName.Data(),fUseReco ? "rec" : "gen")); | |
214 | } | |
215 | ||
ae040693 | 216 | // output |
217 | Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName()); | |
218 | fmyOutput = new TList(); | |
219 | fmyOutput->SetOwner(); | |
220 | fmyOutput->SetName("pippo"); | |
221 | // define histograms | |
222 | DefineHistoForAnalysis(); | |
223 | PostData(1,fmyOutput); | |
224 | ||
225 | fmyOutputF = new TList(); | |
226 | fmyOutputF->SetOwner(); | |
227 | fmyOutputF->SetName("pluto"); | |
228 | // define histograms | |
229 | DefineHistoForFiltering(); | |
230 | PostData(3,fmyOutputF); | |
231 | ||
232 | return; | |
233 | } | |
234 | ||
235 | //_______________________________________________________________________________ | |
236 | ||
237 | void AliAnalysisTaskFlavourFilterAndJetCorrelations::UserExec(Option_t *) | |
238 | { | |
239 | // user exec | |
240 | if (!fInitialized){ | |
241 | AliAnalysisTaskEmcalJet::ExecOnce(); | |
242 | } | |
243 | ||
244 | // Load the event | |
245 | AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent); | |
246 | ||
247 | TClonesArray *arrayDStartoD0pi=0; | |
248 | TClonesArray *trackArr = 0; | |
249 | ||
250 | if (!aodEvent && AODEvent() && IsStandardAOD()) { | |
251 | ||
252 | // In case there is an AOD handler writing a standard AOD, use the AOD | |
253 | // event in memory rather than the input (ESD) event. | |
254 | aodEvent = dynamic_cast<AliAODEvent*> (AODEvent()); | |
255 | ||
256 | // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root) | |
257 | // have to taken from the AOD event hold by the AliAODExtension | |
258 | AliAODHandler* aodHandler = (AliAODHandler*) | |
259 | ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler()); | |
260 | if(aodHandler->GetExtensions()) { | |
261 | AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root"); | |
262 | AliAODEvent *aodFromExt = ext->GetAOD(); | |
263 | arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data()); | |
264 | } | |
265 | } else { | |
266 | arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data()); | |
267 | } | |
268 | ||
269 | if (!arrayDStartoD0pi) { | |
270 | AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data())); | |
271 | // return; | |
272 | } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast())); | |
273 | ||
274 | TClonesArray* mcArray = 0x0; | |
275 | if (fUseMCInfo) { | |
276 | mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName())); | |
277 | if (!mcArray) { | |
278 | printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n"); | |
279 | return; | |
280 | } | |
281 | } | |
282 | ||
283 | //retrieve jets | |
284 | trackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoTracks")); | |
285 | //clusArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClustersCorr")); | |
286 | //jetArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName)); | |
287 | fJetRadius=GetJetContainer(0)->GetJetRadius(); | |
288 | ||
289 | if(!trackArr){ | |
290 | AliInfo(Form("Could not find the track array\n")); | |
291 | return; | |
292 | } | |
293 | ||
294 | ||
295 | //FILTER D MESONS | |
963e733a | 296 | |
297 | fCandidateArray->Clear(); | |
298 | fSideBandArray->Clear(); | |
ae040693 | 299 | |
300 | //Histograms | |
301 | TH1I* hstatF = (TH1I*)fmyOutputF->FindObject("hstatF"); | |
302 | TH1F* hnSBCandEv=(TH1F*)fmyOutputF->FindObject("hnSBCandEv"); | |
303 | TH1F* hnCandEv=(TH1F*)fmyOutputF->FindObject("hnCandEv"); | |
304 | TH2F* hInvMassptDF = (TH2F*)fmyOutputF->FindObject("hInvMassptDF"); | |
305 | ||
306 | TH1F* hPtPion=0x0; | |
307 | if (fCandidateType==kDstartoKpipi) hPtPion = (TH1F*)fmyOutputF->FindObject("hPtPion"); | |
308 | hstatF->Fill(0); | |
309 | ||
310 | // fix for temporary bug in ESDfilter | |
311 | // the AODs with null vertex pointer didn't pass the PhysSel | |
312 | if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return; | |
313 | ||
314 | //Event selection | |
315 | Bool_t iseventselected=fCuts->IsEventSelected(aodEvent); | |
316 | //TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses(); | |
317 | if (!iseventselected) return; | |
318 | hstatF->Fill(1); | |
319 | ||
320 | const Int_t nD = arrayDStartoD0pi->GetEntriesFast(); | |
321 | if(fUseReco) hstatF->Fill(2, nD); | |
322 | ||
323 | //D* and D0 prongs needed to MatchToMC method | |
324 | Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi | |
325 | Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi | |
326 | ||
327 | //D0 from D0 bar | |
328 | Int_t pdgdaughtersD0[2] = { 211, 321 }; // pi,K | |
329 | Int_t pdgdaughtersD0bar[2] = { 321, 211 }; // K,pi | |
330 | ||
331 | Int_t iCand =0; | |
332 | Int_t iSBCand=0; | |
333 | Int_t isSelected = 0; | |
334 | AliAODRecoDecayHF *charmCand = 0; | |
335 | AliAODMCParticle *charmPart = 0; | |
336 | Bool_t isMCBkg=kFALSE; | |
337 | ||
338 | Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass(); | |
339 | ||
340 | Int_t mcLabel = -9999; | |
341 | Int_t pdgMeson = 413; | |
342 | if (fCandidateType==kD0toKpi) pdgMeson = 421; | |
343 | ||
344 | for (Int_t icharm=0; icharm<nD; icharm++) { //loop over D candidates | |
345 | charmCand = (AliAODRecoDecayHF*)arrayDStartoD0pi->At(icharm); // D candidates | |
346 | if (!charmCand) continue; | |
347 | ||
348 | ||
349 | if (fUseMCInfo) { // Look in MC, try to simulate the z | |
350 | if (fCandidateType==kDstartoKpipi) { | |
351 | AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand; | |
352 | mcLabel = temp->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray); | |
353 | } | |
354 | ||
355 | if (fCandidateType==kD0toKpi) | |
356 | mcLabel = charmCand->MatchToMC(421,mcArray,fNProngs,fPDGdaughters); | |
357 | ||
358 | if (mcLabel<=0) isMCBkg=kTRUE; | |
359 | else hstatF->Fill(2); | |
360 | if (!isMCBkg) charmPart=(AliAODMCParticle*)mcArray->At(mcLabel); | |
361 | } | |
362 | ||
363 | Double_t ptD = charmCand->Pt(); | |
364 | ||
365 | // region of interest + cuts | |
366 | if (!fCuts->IsInFiducialAcceptance(ptD,charmCand->Y(pdgMeson))) continue; | |
367 | ||
368 | if(!fUseMCInfo && fCandidateType==kDstartoKpipi){ | |
369 | //select by track cuts the side band candidates (don't want mass cut) | |
370 | isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kTracks,aodEvent); | |
371 | if (!isSelected) continue; | |
372 | //add a reasonable cut on the invariant mass (e.g. (+-2\sigma, +-10 \sigma), with \sigma = fSigmaD0[bin]) | |
373 | Int_t bin = fCuts->PtBin(ptD); | |
374 | if(bin<0 || bin>=fCuts->GetNPtBins()) { | |
375 | AliError(Form("Pt %.3f out of bounds",ptD)); | |
376 | continue; | |
377 | } | |
378 | AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand; | |
379 | //if data and Dstar from D0 side band | |
380 | if (((temp->InvMassD0()<=(mPDGD0-3.*fSigmaD0[bin])) && (temp->InvMassD0()>(mPDGD0-10.*fSigmaD0[bin]))) /*left side band*/|| ((temp->InvMassD0()>=(mPDGD0+3.*fSigmaD0[bin])) && (temp->InvMassD0()<(mPDGD0+10.*fSigmaD0[bin])))/*right side band*/){ | |
381 | ||
382 | new ((*fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(*temp); | |
383 | iSBCand++; | |
384 | } | |
385 | } | |
386 | //candidate selected by cuts and PID | |
387 | isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kAll,aodEvent); //selected | |
388 | if (!isSelected) continue; | |
389 | ||
390 | //for data and MC signal fill fCandidateArray | |
391 | if(!fUseMCInfo || (fUseMCInfo && !isMCBkg)){ | |
392 | // for data or MC with the requirement fUseReco fill with candidates | |
393 | if(fUseReco) { | |
394 | new ((*fCandidateArray)[iCand]) AliAODRecoDecayHF(*charmCand); | |
395 | //Printf("Filling reco"); | |
396 | hstatF->Fill(3); | |
397 | } | |
398 | // for MC with requirement particle level fill with AliAODMCParticle | |
399 | else if (fUseMCInfo) { | |
400 | new ((*fCandidateArray)[iCand]) AliAODMCParticle(*charmPart); | |
401 | //Printf("Filling gen"); | |
402 | hstatF->Fill(3); | |
403 | } | |
404 | ||
405 | iCand++; | |
406 | } | |
407 | //for MC background fill fSideBandArray (which is instead filled above for DStar in case of data for the side bands candidates) | |
408 | else if(fUseReco){ | |
409 | new ((*fSideBandArray)[iSBCand]) AliAODRecoDecayHF(*charmCand); | |
410 | iSBCand++; | |
411 | } | |
412 | ||
413 | ||
414 | Double_t masses[2]; | |
415 | if (fCandidateType==kDstartoKpipi) { //D*->D0pi->Kpipi | |
416 | ||
417 | //softpion from D* decay | |
418 | AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand; | |
419 | AliAODTrack *track2 = (AliAODTrack*)temp->GetBachelor(); | |
420 | ||
421 | // select D* in the D0 window. | |
422 | // In the cut object window is loose to allow for side bands | |
423 | ||
424 | ||
425 | // retrieve the corresponding pt bin for the candidate | |
426 | // and set the expected D0 width (x3) | |
427 | // static const Int_t n = fCuts->GetNPtBins(); | |
428 | Int_t bin = fCuts->PtBin(ptD); | |
429 | if(bin<0 || bin>=fCuts->GetNPtBins()) { | |
430 | AliError(Form("Pt %.3f out of bounds",ptD)); | |
431 | continue; | |
432 | } | |
433 | ||
434 | AliInfo(Form("Pt bin %d and sigma D0 %.4f",bin,fSigmaD0[bin])); | |
435 | //consider the Dstar candidates only if the mass of the D0 is in 3 sigma wrt the PDG value | |
436 | if ((temp->InvMassD0()>=(mPDGD0-3.*fSigmaD0[bin])) && (temp->InvMassD0()<=(mPDGD0+3.*fSigmaD0[bin]))) { | |
437 | masses[0] = temp->DeltaInvMass(); //D* | |
438 | masses[1] = 0.; //dummy for D* | |
439 | ||
440 | //D* delta mass | |
441 | hInvMassptDF->Fill(masses[0], ptD); // 2 D slice for pt bins | |
442 | ||
443 | // D* pt and soft pion pt for good candidates | |
444 | Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass(); | |
445 | Double_t invmassDelta = temp->DeltaInvMass(); | |
446 | if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0021) hPtPion->Fill(track2->Pt()); | |
447 | } | |
448 | } //Dstar specific | |
449 | ||
450 | if (fCandidateType==kD0toKpi) { //D0->Kpi | |
451 | ||
452 | //needed quantities | |
453 | masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0 | |
454 | masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar | |
455 | hstatF->Fill(3); | |
456 | ||
457 | // mass vs pt | |
458 | if (isSelected==1 || isSelected==3) hInvMassptDF->Fill(masses[0],ptD); | |
459 | if (isSelected>=2) hInvMassptDF->Fill(masses[1],ptD); | |
460 | } //D0 specific | |
461 | ||
462 | charmCand = 0; | |
463 | isMCBkg=kFALSE; | |
464 | } // end of D cand loop | |
465 | ||
466 | hnCandEv->Fill(fCandidateArray->GetEntriesFast()); | |
467 | if (fCandidateType==kDstartoKpipi || fUseMCInfo) { | |
468 | Int_t nsbcand=fSideBandArray->GetEntriesFast(); | |
469 | hstatF->Fill(4,nsbcand); | |
470 | hnSBCandEv->Fill(nsbcand); | |
471 | } | |
472 | ||
473 | //CORRELATION WITH JETS | |
474 | /* | |
475 | TString arrname="fCandidateArray"; | |
476 | TString candarrname=Form("%s%s%s",arrname.Data(),fCandArrName.Data(),fUseReco ? "rec" : "gen"); | |
477 | fCandidateArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(candarrname)); | |
478 | if (!fCandidateArray) { | |
479 | Printf("%s array not found",candarrname.Data()); | |
480 | InputEvent()->GetList()->ls(); | |
481 | return; | |
482 | } | |
483 | //Printf("ncandidates found %d",fCandidateArray->GetEntriesFast()); | |
484 | ||
485 | TString arrSBname="fSideBandArray"; | |
486 | TString sbarrname=Form("%s%s%s",arrSBname.Data(),fCandArrName.Data(),fUseReco ? "rec" : "gen"); | |
487 | fSideBandArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sbarrname)); | |
488 | if (fCandidateType==1 && !fSideBandArray) { | |
489 | Printf("%s array not found",sbarrname.Data()); | |
490 | InputEvent()->GetList()->ls(); | |
491 | return; | |
492 | } | |
493 | //Printf("nSBCand or Bkg found %d",fSideBandArray->GetEntriesFast()); | |
494 | */ | |
495 | ||
496 | //Histograms | |
497 | TH1I* hstat=(TH1I*)fmyOutput->FindObject("hstat"); | |
498 | TH1F* hPtJetTrks=(TH1F*)fmyOutput->FindObject("hPtJetTrks"); | |
499 | TH1F* hPhiJetTrks=(TH1F*)fmyOutput->FindObject("hPhiJetTrks"); | |
500 | TH1F* hEtaJetTrks=(TH1F*)fmyOutput->FindObject("hEtaJetTrks"); | |
501 | TH1F* hEjetTrks=(TH1F*)fmyOutput->FindObject("hEjetTrks"); | |
502 | TH1F* hPtJet=(TH1F*)fmyOutput->FindObject("hPtJet"); | |
503 | TH1F* hPhiJet=(TH1F*)fmyOutput->FindObject("hPhiJet"); | |
504 | TH1F* hEtaJet=(TH1F*)fmyOutput->FindObject("hEtaJet"); | |
505 | TH1F* hEjet=(TH1F*)fmyOutput->FindObject("hEjet"); | |
506 | TH1F* hNtrArr=(TH1F*)fmyOutput->FindObject("hNtrArr"); | |
507 | TH1F* hNJetPerEv=(TH1F*)fmyOutput->FindObject("hNJetPerEv"); | |
508 | TH1F* hdeltaRJetTracks=(TH1F*)fmyOutput->FindObject("hdeltaRJetTracks"); | |
509 | TH1F* hNDPerEvNoJet=(TH1F*)fmyOutput->FindObject("hNDPerEvNoJet"); | |
510 | TH1F* hptDPerEvNoJet=(TH1F*)fmyOutput->FindObject("hptDPerEvNoJet"); | |
511 | TH1F* hNJetPerEvNoD=(TH1F*)fmyOutput->FindObject("hNJetPerEvNoD"); | |
512 | TH1F* hPtJetPerEvNoD=(TH1F*)fmyOutput->FindObject("hPtJetPerEvNoD"); | |
513 | THnSparseF* hnspDstandalone=(THnSparseF*)fmyOutput->FindObject("hsDstandalone"); | |
514 | ||
515 | hstat->Fill(0); | |
516 | ||
517 | // fix for temporary bug in ESDfilter | |
518 | // the AODs with null vertex pointer didn't pass the PhysSel | |
519 | if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return; | |
520 | ||
521 | //Event selection | |
522 | //Bool_t iseventselected=fCuts->IsEventSelected(aodEvent); | |
523 | TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses(); | |
524 | if(!iseventselected) return; | |
525 | ||
526 | hstat->Fill(1); | |
527 | ||
528 | //retrieve charm candidates selected | |
529 | Int_t candidates = fCandidateArray->GetEntriesFast(); | |
530 | ||
531 | //trigger on jets | |
532 | ||
533 | Int_t njets=GetJetContainer()->GetNJets(); | |
534 | if(njets == 0) { | |
535 | hstat->Fill(6, candidates); | |
536 | hNDPerEvNoJet->Fill(candidates); | |
537 | for(Int_t iD=0;iD<candidates;iD++){ | |
538 | AliVParticle* cand=(AliVParticle*)fCandidateArray->At(iD); | |
539 | if(!cand) continue; | |
540 | hptDPerEvNoJet->Fill(cand->Pt()); | |
541 | ||
542 | } | |
543 | return; | |
544 | ||
545 | } | |
546 | ||
547 | //loop on candidates standalone (checking the candidates are there and their phi-eta distributions) | |
548 | ||
549 | for(Int_t ic = 0; ic < candidates; ic++) { | |
550 | ||
551 | // D* candidates | |
552 | AliVParticle* charm=0x0; | |
553 | charm=(AliVParticle*)fCandidateArray->At(ic); | |
554 | if(!charm) continue; | |
555 | hstat->Fill(2); | |
556 | ||
557 | Double_t candsparse[4]={charm->Eta(), charm->Phi(), charm->Pt(), 0}; | |
558 | ||
559 | if(fCandidateType==kDstartoKpipi) { | |
560 | AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)charm; | |
561 | Double_t deltamass= dstar->DeltaInvMass(); | |
562 | candsparse[3]=deltamass; | |
563 | hnspDstandalone->Fill(candsparse); | |
564 | } | |
565 | if(fCandidateType==kD0toKpi){ | |
566 | AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)charm; | |
567 | Int_t isselected=fCuts->IsSelected(dzero,AliRDHFCuts::kAll,aodEvent); | |
568 | ||
569 | Double_t masses[2]; | |
570 | //Int_t pdgdaughtersD0[2]={211,321};//pi,K | |
571 | //Int_t pdgdaughtersD0bar[2]={321,211};//K,pi | |
572 | ||
573 | masses[0]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0 | |
574 | masses[1]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar | |
575 | if(isselected==1 || isselected==3) { | |
576 | candsparse[3]=masses[0]; | |
577 | hnspDstandalone->Fill(candsparse); | |
578 | } | |
579 | if(isselected>=2){ | |
580 | candsparse[3]=masses[1]; | |
581 | hnspDstandalone->Fill(candsparse); | |
582 | ||
583 | } | |
584 | } | |
585 | } | |
586 | ||
587 | // we start with jets | |
588 | Double_t ejet = 0; | |
589 | Double_t phiJet = 0; | |
590 | Double_t etaJet = 0; | |
591 | Double_t ptjet = 0; | |
592 | Double_t leadingJet =0; | |
593 | ||
594 | Int_t ntrarr=trackArr->GetEntriesFast(); | |
595 | hNtrArr->Fill(ntrarr); | |
596 | ||
597 | for(Int_t i=0;i<ntrarr;i++){ | |
598 | AliVTrack *jtrack=static_cast<AliVTrack*>(trackArr->At(i)); | |
599 | //reusing histograms | |
600 | hPtJetTrks->Fill(jtrack->Pt()); | |
601 | hPhiJetTrks->Fill(jtrack->Phi()); | |
602 | hEtaJetTrks->Fill(jtrack->Eta()); | |
603 | hEjetTrks->Fill(jtrack->E()); | |
604 | } | |
605 | ||
606 | ||
607 | //option to use only the leading jet | |
608 | if(fLeadingJetOnly){ | |
609 | for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) { | |
610 | AliEmcalJet* jetL = (AliEmcalJet*)GetJetFromArray(iJetsL); | |
611 | ptjet = jetL->Pt(); | |
612 | if(ptjet>leadingJet ) leadingJet = ptjet; | |
613 | ||
614 | } | |
615 | } | |
616 | ||
617 | Int_t cntjet=0; | |
618 | //loop over jets and consider only the leading jet in the event | |
619 | for (Int_t iJets = 0; iJets<njets; iJets++) { | |
620 | //Printf("Jet N %d",iJets); | |
621 | AliEmcalJet* jet = (AliEmcalJet*)GetJetFromArray(iJets); | |
622 | //Printf("Corr task Accept Jet"); | |
623 | if(!AcceptJet(jet)) { | |
624 | hstat->Fill(5); | |
625 | continue; | |
626 | } | |
627 | ||
628 | //jets variables | |
629 | ejet = jet->E(); | |
630 | phiJet = jet->Phi(); | |
631 | etaJet = jet->Eta(); | |
632 | ptjet = jet->Pt(); | |
633 | ||
634 | // choose the leading jet | |
635 | if(fLeadingJetOnly && (ejet<leadingJet)) continue; | |
636 | //used for normalization | |
637 | hstat->Fill(3); | |
638 | cntjet++; | |
639 | // fill energy, eta and phi of the jet | |
640 | hEjet ->Fill(ejet); | |
641 | hPhiJet ->Fill(phiJet); | |
642 | hEtaJet ->Fill(etaJet); | |
643 | hPtJet ->Fill(ptjet); | |
644 | ||
645 | //loop on jet particles | |
646 | Int_t ntrjet= jet->GetNumberOfTracks(); | |
647 | for(Int_t itrk=0;itrk<ntrjet;itrk++){ | |
648 | ||
649 | AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,trackArr); | |
650 | hdeltaRJetTracks->Fill(DeltaR(jet,jetTrk)); | |
651 | ||
652 | }//end loop on jet tracks | |
653 | ||
654 | if(candidates==0){ | |
655 | hstat->Fill(7); | |
656 | hPtJetPerEvNoD->Fill(jet->Pt()); | |
657 | } | |
658 | ||
659 | //Printf("N candidates %d ", candidates); | |
660 | for(Int_t ic = 0; ic < candidates; ic++) { | |
661 | ||
662 | // D* candidates | |
663 | AliVParticle* charm=0x0; | |
664 | charm=(AliVParticle*)fCandidateArray->At(ic); | |
665 | if(!charm) continue; | |
666 | ||
667 | ||
668 | FlagFlavour(charm, jet); | |
669 | if (jet->TestFlavourTag(AliEmcalJet::kDStar)) hstat->Fill(4); | |
670 | ||
671 | FillHistogramsRecoJetCorr(charm, jet); | |
672 | ||
673 | } | |
674 | //retrieve side band background candidates for Dstar | |
675 | Int_t nsbcand = 0; if(fSideBandArray) nsbcand=fSideBandArray->GetEntriesFast(); | |
676 | ||
677 | for(Int_t ib=0;ib<nsbcand;ib++){ | |
678 | if(fCandidateType==kDstartoKpipi && !fUseMCInfo){ | |
679 | AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)fSideBandArray->At(ib); | |
680 | if(!sbcand) continue; | |
681 | SideBandBackground(sbcand,jet); | |
682 | } | |
683 | if(fUseMCInfo){ | |
684 | AliAODRecoDecayHF* charmbg = 0x0; | |
685 | charmbg=(AliAODRecoDecayHF*)fCandidateArray->At(ib); | |
686 | if(!charmbg) continue; | |
687 | MCBackground(charmbg,jet); | |
688 | } | |
689 | } | |
690 | } // end of jet loop | |
691 | ||
692 | hNJetPerEv->Fill(cntjet); | |
693 | if(candidates==0) hNJetPerEvNoD->Fill(cntjet); | |
694 | PostData(1,fmyOutput); | |
695 | PostData(3,fmyOutputF); | |
696 | } | |
697 | ||
698 | //_______________________________________________________________________________ | |
699 | ||
700 | void AliAnalysisTaskFlavourFilterAndJetCorrelations::Terminate(Option_t*) | |
701 | { | |
702 | // The Terminate() function is the last function to be called during | |
703 | // a query. It always runs on the client, it can be used to present | |
704 | // the results graphically or save the results to file. | |
705 | ||
706 | Info("Terminate"," terminate"); | |
707 | AliAnalysisTaskSE::Terminate(); | |
708 | ||
709 | fmyOutput = dynamic_cast<TList*> (GetOutputData(1)); | |
710 | if (!fmyOutput) { | |
711 | printf("ERROR: fmyOutput not available\n"); | |
712 | return; | |
713 | } | |
714 | } | |
715 | ||
716 | //_______________________________________________________________________________ | |
717 | ||
718 | void AliAnalysisTaskFlavourFilterAndJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){ | |
719 | Float_t mass=0; | |
720 | Int_t abspdg=TMath::Abs(pdg); | |
721 | ||
722 | mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass(); | |
723 | // compute the Delta mass for the D* | |
724 | if(fCandidateType==kDstartoKpipi){ | |
725 | Float_t mass1=0; | |
726 | mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass(); | |
727 | mass = mass-mass1; | |
728 | } | |
729 | ||
730 | fMinMass = mass-range; | |
731 | fMaxMass = mass+range; | |
732 | ||
733 | AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass)); | |
734 | if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n"); | |
735 | } | |
736 | ||
737 | //_______________________________________________________________________________ | |
738 | ||
739 | void AliAnalysisTaskFlavourFilterAndJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){ | |
740 | if(uplimit>lowlimit) | |
741 | { | |
742 | fMinMass = lowlimit; | |
743 | fMaxMass = uplimit; | |
744 | } | |
745 | else{ | |
746 | printf("Error! Lower limit larger than upper limit!\n"); | |
747 | fMinMass = uplimit - uplimit*0.2; | |
748 | fMaxMass = uplimit; | |
749 | } | |
750 | } | |
751 | ||
752 | //_______________________________________________________________________________ | |
753 | ||
754 | Bool_t AliAnalysisTaskFlavourFilterAndJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){ | |
755 | if(nptbins>30) { | |
756 | AliInfo("Maximum number of bins allowed is 30!"); | |
757 | return kFALSE; | |
758 | } | |
759 | if(!width) return kFALSE; | |
760 | for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt]; | |
761 | return kTRUE; | |
762 | } | |
763 | ||
764 | //_______________________________________________________________________________ | |
765 | ||
766 | Double_t AliAnalysisTaskFlavourFilterAndJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet) const{ | |
767 | if(!part || !jet){ | |
768 | printf("Particle or jet do not exist!\n"); | |
769 | return -99; | |
770 | } | |
771 | Double_t p[3],pj[3]; | |
772 | Bool_t okpp=part->PxPyPz(p); | |
773 | Bool_t okpj=jet->PxPyPz(pj); | |
774 | if(!okpp || !okpj){ | |
775 | printf("Problems getting momenta\n"); | |
776 | return -999; | |
777 | } | |
778 | Double_t pjet=jet->P(); | |
779 | Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet*pjet); | |
780 | return z; | |
781 | } | |
782 | ||
783 | //_______________________________________________________________________________ | |
784 | ||
785 | Bool_t AliAnalysisTaskFlavourFilterAndJetCorrelations::DefineHistoForAnalysis(){ | |
786 | ||
787 | // Statistics | |
788 | TH1I* hstat=new TH1I("hstat","Statistics",8,-0.5,7.5); | |
789 | hstat->GetXaxis()->SetBinLabel(1,"N ev anal"); | |
790 | hstat->GetXaxis()->SetBinLabel(2,"N ev sel"); | |
791 | hstat->GetXaxis()->SetBinLabel(3,"N cand sel & jet"); | |
792 | hstat->GetXaxis()->SetBinLabel(4,"N jets"); | |
793 | hstat->GetXaxis()->SetBinLabel(5,"N cand in jet"); | |
794 | hstat->GetXaxis()->SetBinLabel(6,"N jet rej"); | |
795 | hstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet"); | |
796 | hstat->GetXaxis()->SetBinLabel(8,"N jets & !D"); | |
797 | hstat->SetNdivisions(1); | |
798 | fmyOutput->Add(hstat); | |
799 | ||
800 | const Int_t nbinsmass=200; | |
801 | const Int_t nbinsptjet=500; | |
802 | const Int_t nbinsptD=100; | |
803 | const Int_t nbinsz=100; | |
804 | const Int_t nbinsphi=200; | |
805 | const Int_t nbinseta=100; | |
806 | ||
807 | const Float_t ptjetlims[2]={0.,200.}; | |
808 | const Float_t ptDlims[2]={0.,50.}; | |
809 | const Float_t zlims[2]={0.,1.2}; | |
810 | const Float_t philims[2]={0.,6.3}; | |
811 | const Float_t etalims[2]={-1.5,1.5}; | |
812 | ||
813 | if(fCandidateType==kDstartoKpipi) | |
814 | { | |
815 | ||
816 | TH2F* hDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]); | |
817 | hDiffSideBand->SetStats(kTRUE); | |
818 | hDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV"); | |
819 | hDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)"); | |
820 | hDiffSideBand->Sumw2(); | |
821 | fmyOutput->Add(hDiffSideBand); | |
822 | ||
823 | ||
824 | TH1F* hPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10); | |
825 | hPtPion->SetStats(kTRUE); | |
826 | hPtPion->GetXaxis()->SetTitle("GeV/c"); | |
827 | hPtPion->GetYaxis()->SetTitle("Entries"); | |
828 | hPtPion->Sumw2(); | |
829 | fmyOutput->Add(hPtPion); | |
830 | ||
831 | } | |
832 | // D related histograms | |
833 | TH1F *hNDPerEvNoJet=new TH1F("hNDPerEvNoJet","Number of candidates per event with no jets; N candidate/ev with no jet", 20, 0., 20.); | |
834 | hNDPerEvNoJet->Sumw2(); | |
835 | fmyOutput->Add(hNDPerEvNoJet); | |
836 | ||
837 | TH1F *hptDPerEvNoJet=new TH1F("hptDPerEvNoJet","pt distribution of candidates per events with no jets; p_{t}^{D} (GeV/c)",nbinsptD, ptDlims[0],ptDlims[1]); | |
838 | hptDPerEvNoJet->Sumw2(); | |
839 | fmyOutput->Add(hptDPerEvNoJet); | |
840 | ||
841 | const Int_t nAxisD=4; | |
842 | const Int_t nbinsSparseD[nAxisD]={nbinsphi,nbinseta,nbinsptD,nbinsmass}; | |
843 | const Double_t minSparseD[nAxisD] ={philims[0],etalims[0],ptDlims[0],fMinMass}; | |
844 | const Double_t maxSparseD[nAxisD] ={philims[1],etalims[1],ptDlims[1],fMaxMass}; | |
845 | THnSparseF *hsDstandalone=new THnSparseF("hsDstandalone","#phi, #eta, p_{T}^{D}, and mass", nAxisD, nbinsSparseD, minSparseD, maxSparseD); | |
846 | hsDstandalone->Sumw2(); | |
847 | ||
848 | fmyOutput->Add(hsDstandalone); | |
849 | ||
850 | // jet related fistograms | |
851 | ||
852 | TH1F* hEjetTrks = new TH1F("hEjetTrks", "Jet tracks energy distribution;Energy (GeV)",500,0,200); | |
853 | hEjetTrks->Sumw2(); | |
854 | TH1F* hPhiJetTrks = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]); | |
855 | hPhiJetTrks->Sumw2(); | |
856 | TH1F* hEtaJetTrks = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta", nbinseta,etalims[0],etalims[1]); | |
857 | hEtaJetTrks->Sumw2(); | |
858 | TH1F* hPtJetTrks = new TH1F("hPtJetTrks", "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]); | |
859 | hPtJetTrks->Sumw2(); | |
860 | ||
861 | TH1F* hEjet = new TH1F("hEjet", "Jet energy distribution;Energy (GeV)",500,0,200); | |
862 | hEjet->Sumw2(); | |
863 | TH1F* hPhiJet = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]); | |
864 | hPhiJet->Sumw2(); | |
865 | TH1F* hEtaJet = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]); | |
866 | hEtaJet->Sumw2(); | |
867 | TH1F* hPtJet = new TH1F("hPtJet", "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]); | |
868 | hPtJet->Sumw2(); | |
869 | ||
870 | 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]); | |
871 | hPtJetWithD->Sumw2(); | |
872 | //for the MC this histogram is filled with the real background | |
873 | 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]); | |
874 | hPtJetWithDsb->Sumw2(); | |
875 | ||
876 | TH1F* hdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.); | |
877 | hdeltaRJetTracks->Sumw2(); | |
878 | ||
879 | TH1F* hNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000); | |
880 | hNtrArr->Sumw2(); | |
881 | ||
882 | TH1F *hNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5); | |
883 | hNJetPerEv->Sumw2(); | |
884 | ||
885 | 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); | |
886 | hNJetPerEvNoD->Sumw2(); | |
887 | ||
888 | TH1F *hPtJetPerEvNoD=new TH1F("hPtJetPerEvNoD","pt distribution of jets per event with no D; p_{T}^{jet} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]); | |
889 | hPtJetPerEvNoD->Sumw2(); | |
890 | ||
891 | fmyOutput->Add(hEjetTrks); | |
892 | fmyOutput->Add(hPhiJetTrks); | |
893 | fmyOutput->Add(hEtaJetTrks); | |
894 | fmyOutput->Add(hPtJetTrks); | |
895 | fmyOutput->Add(hEjet); | |
896 | fmyOutput->Add(hPhiJet); | |
897 | fmyOutput->Add(hEtaJet); | |
898 | fmyOutput->Add(hPtJet); | |
899 | fmyOutput->Add(hPtJetWithD); | |
900 | fmyOutput->Add(hPtJetWithDsb); | |
901 | fmyOutput->Add(hdeltaRJetTracks); | |
902 | fmyOutput->Add(hNtrArr); | |
903 | fmyOutput->Add(hNJetPerEv); | |
904 | fmyOutput->Add(hNJetPerEvNoD); | |
905 | fmyOutput->Add(hPtJetPerEvNoD); | |
906 | ||
907 | TH1F* hDeltaRD=new TH1F("hDeltaRD","#Delta R distribution of D candidates selected;#Delta R",200, 0.,10.); | |
908 | hDeltaRD->Sumw2(); | |
909 | fmyOutput->Add(hDeltaRD); | |
910 | ||
911 | //background (side bands for the Dstar and like sign for D0) | |
912 | fJetRadius=GetJetContainer(0)->GetJetRadius(); | |
913 | 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]); | |
914 | hInvMassptD->SetStats(kTRUE); | |
915 | hInvMassptD->GetXaxis()->SetTitle("mass (GeV)"); | |
916 | hInvMassptD->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)"); | |
917 | hInvMassptD->Sumw2(); | |
918 | ||
919 | fmyOutput->Add(hInvMassptD); | |
920 | ||
921 | if(fUseMCInfo){ | |
922 | 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]); | |
923 | hInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)")); | |
924 | hInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)"); | |
925 | hInvMassptDbg->Sumw2(); | |
926 | fmyOutput->Add(hInvMassptDbg); | |
927 | ||
928 | } | |
929 | Double_t pi=TMath::Pi(), philims2[2]; | |
930 | philims2[0]=-pi*1./2.; | |
931 | philims2[1]=pi*3./2.; | |
932 | const Int_t nAxis=6; | |
933 | const Int_t nbinsSparse[nAxis]={nbinsz,nbinsphi,nbinsptjet,nbinsptD,nbinsmass,2}; | |
934 | const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5}; | |
935 | const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5}; | |
936 | THnSparseF *hsDphiz=new THnSparseF("hsDphiz","Z and #Delta#phi vs p_{T}^{jet}, p_{T}^{D}, and mass", nAxis, nbinsSparse, minSparse, maxSparse); | |
937 | hsDphiz->Sumw2(); | |
938 | ||
939 | fmyOutput->Add(hsDphiz); | |
940 | ||
941 | PostData(1, fmyOutput); | |
942 | ||
943 | return kTRUE; | |
944 | } | |
945 | ||
946 | //_______________________________________________________________________________ | |
947 | ||
948 | Bool_t AliAnalysisTaskFlavourFilterAndJetCorrelations::DefineHistoForFiltering() | |
949 | { | |
950 | // | |
951 | // AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis | |
952 | // | |
953 | ||
954 | // Statistics | |
955 | TH1I* hstat = new TH1I("hstatF","Statistics",5,-0.5,4.5); | |
956 | hstat->GetXaxis()->SetBinLabel(1, "N ev anal"); | |
957 | hstat->GetXaxis()->SetBinLabel(2, "N ev sel"); | |
958 | if(fUseReco) hstat->GetXaxis()->SetBinLabel(3, "N cand"); | |
959 | else hstat->GetXaxis()->SetBinLabel(3, "N Gen D"); | |
960 | hstat->GetXaxis()->SetBinLabel(4, "N cand sel cuts"); | |
961 | if (fCandidateType==kDstartoKpipi) hstat->GetXaxis()->SetBinLabel(5, "N side band cand"); | |
962 | hstat->SetNdivisions(1); | |
963 | fmyOutputF->Add(hstat); | |
964 | ||
965 | TH1F* hnCandEv=new TH1F("hnCandEv", "Number of candidates per event (after cuts);# cand/ev", 100, 0.,100.); | |
966 | fmyOutputF->Add(hnCandEv); | |
967 | ||
968 | // Invariant mass related histograms | |
969 | const Int_t nbinsmass = 200; | |
970 | TH2F *hInvMass = new TH2F("hInvMassptDF", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, 100, 0., 50.); | |
971 | hInvMass->SetStats(kTRUE); | |
972 | hInvMass->GetXaxis()->SetTitle("mass (GeV/c)"); | |
973 | hInvMass->GetYaxis()->SetTitle("p_{T} (GeV/c)"); | |
974 | fmyOutputF->Add(hInvMass); | |
975 | ||
976 | if (fCandidateType==kDstartoKpipi) { | |
977 | TH1F* hnSBCandEv=new TH1F("hnSBCandEv", "Number of side bands candidates per event (after cuts);# cand/ev", 100, 0.,100.); | |
978 | fmyOutputF->Add(hnSBCandEv); | |
979 | ||
980 | TH1F* hPtPion = new TH1F("hPtPion", "Primary pions candidates pt", 500, 0., 10.); | |
981 | hPtPion->SetStats(kTRUE); | |
982 | hPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)"); | |
983 | hPtPion->GetYaxis()->SetTitle("entries"); | |
984 | fmyOutputF->Add(hPtPion); | |
985 | } | |
986 | PostData(3, fmyOutputF); | |
987 | return kTRUE; | |
988 | } | |
989 | ||
990 | //_______________________________________________________________________________ | |
991 | ||
992 | void AliAnalysisTaskFlavourFilterAndJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet){ | |
993 | ||
994 | Double_t ptD=candidate->Pt(); | |
995 | Double_t ptjet=jet->Pt(); | |
996 | Double_t deltaR=DeltaR(candidate,jet); | |
997 | Double_t deltaphi = jet->Phi()-candidate->Phi(); | |
998 | if(deltaphi<=-(TMath::Pi())/2) deltaphi = deltaphi+2*(TMath::Pi()); | |
999 | if(deltaphi>(3*(TMath::Pi()))/2) deltaphi = deltaphi-2*(TMath::Pi()); | |
1000 | Double_t z=Z(candidate,jet); | |
1001 | ||
1002 | TH1F* hDeltaRD=(TH1F*)fmyOutput->FindObject("hDeltaRD"); | |
1003 | hDeltaRD->Fill(deltaR); | |
1004 | if(fUseReco){ | |
1005 | if(fCandidateType==kD0toKpi) { | |
1006 | AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate; | |
1007 | FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,ptjet,deltaR, AODEvent()); | |
1008 | ||
1009 | } | |
1010 | ||
1011 | if(fCandidateType==kDstartoKpipi) { | |
1012 | AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate; | |
1013 | FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,ptjet,deltaR); | |
1014 | ||
1015 | } | |
1016 | } else{ //generated level | |
1017 | //AliInfo("Non reco"); | |
1018 | FillHistogramsMCGenDJetCorr(deltaphi,z,ptD,ptjet,deltaR); | |
1019 | } | |
1020 | } | |
1021 | ||
1022 | //_______________________________________________________________________________ | |
1023 | ||
1024 | void AliAnalysisTaskFlavourFilterAndJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj,Double_t deltaR, AliAODEvent* aodEvent){ | |
1025 | ||
1026 | ||
1027 | Double_t masses[2]={0.,0.}; | |
1028 | Int_t pdgdaughtersD0[2]={211,321};//pi,K | |
1029 | Int_t pdgdaughtersD0bar[2]={321,211};//K,pi | |
1030 | ||
1031 | masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0 | |
1032 | masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar | |
1033 | ||
1034 | TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD"); | |
1035 | THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz"); | |
1036 | Double_t point[5]={z,dPhi,ptj,ptD,masses[0]}; | |
1037 | ||
1038 | Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent); | |
1039 | if(isselected==1 || isselected==3) { | |
1040 | ||
1041 | if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[0],ptD); | |
1042 | ||
1043 | FillMassHistograms(masses[0], ptD, deltaR); | |
1044 | hsDphiz->Fill(point,1.); | |
1045 | } | |
1046 | if(isselected>=2) { | |
1047 | if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[1],ptD); | |
1048 | ||
1049 | FillMassHistograms(masses[1], ptD, deltaR); | |
1050 | point[4]=masses[1]; | |
1051 | hsDphiz->Fill(point,1.); | |
1052 | } | |
1053 | ||
1054 | } | |
1055 | ||
1056 | //_______________________________________________________________________________ | |
1057 | ||
1058 | void AliAnalysisTaskFlavourFilterAndJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Double_t deltaR){ | |
1059 | //dPhi and z not used at the moment,but will be (re)added | |
1060 | ||
1061 | AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor(); | |
1062 | Double_t deltamass= dstar->DeltaInvMass(); | |
1063 | //Double_t massD0= dstar->InvMassD0(); | |
1064 | ||
1065 | ||
1066 | TH1F* hPtPion=(TH1F*)fmyOutput->FindObject("hPtPion"); | |
1067 | hPtPion->Fill(softpi->Pt()); | |
1068 | ||
1069 | TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD"); | |
1070 | THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz"); | |
1071 | Int_t isSB=IsDzeroSideBand(dstar); | |
1072 | ||
1073 | Double_t point[6]={z,dPhi,ptj,ptD,deltamass,isSB}; | |
1074 | ||
1075 | if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,deltamass,ptD); | |
1076 | ||
1077 | FillMassHistograms(deltamass, ptD, deltaR); | |
1078 | hsDphiz->Fill(point,1.); | |
1079 | ||
1080 | } | |
1081 | ||
1082 | //_______________________________________________________________________________ | |
1083 | ||
1084 | void AliAnalysisTaskFlavourFilterAndJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t z,Double_t ptD,Double_t ptjet,Double_t deltaR){ | |
1085 | ||
1086 | Double_t pdgmass=0; | |
1087 | TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD"); | |
1088 | THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz"); | |
1089 | Double_t point[6]={z,dPhi,ptjet,ptD,pdgmass,0}; | |
1090 | ||
1091 | if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass(); | |
1092 | if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass(); | |
1093 | point[4]=pdgmass; | |
1094 | ||
1095 | if(deltaR<fJetRadius) { | |
1096 | hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet | |
1097 | hsDphiz->Fill(point,1.); | |
1098 | } | |
1099 | ||
1100 | } | |
1101 | ||
1102 | //_______________________________________________________________________________ | |
1103 | ||
1104 | void AliAnalysisTaskFlavourFilterAndJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD, Double_t deltaR){ | |
1105 | ||
1106 | if(deltaR<fJetRadius) { | |
1107 | TH2F* hInvMassptD=(TH2F*)fmyOutput->FindObject("hInvMassptD"); | |
1108 | hInvMassptD->Fill(mass,ptD); | |
1109 | } | |
1110 | } | |
1111 | ||
1112 | //________________________________________________________________________________ | |
1113 | ||
1114 | void AliAnalysisTaskFlavourFilterAndJetCorrelations::FlagFlavour(AliVParticle *charm, AliEmcalJet *jet){ | |
1115 | Double_t deltaR=DeltaR(charm, jet); | |
1116 | AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar; | |
1117 | if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0; | |
1118 | if (deltaR<fJetRadius) jet->AddFlavourTag(tag); | |
1119 | ||
1120 | return; | |
1121 | ||
1122 | } | |
1123 | ||
1124 | //_______________________________________________________________________________ | |
1125 | ||
1126 | void AliAnalysisTaskFlavourFilterAndJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){ | |
1127 | ||
1128 | // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas | |
1129 | // (expected detector resolution) on the left and right frm the D0 mass. Each band | |
1130 | // has a width of ~5 sigmas. Two band needed for opening angle considerations | |
1131 | TH2F* hDiffSideBand=(TH2F*)fmyOutput->FindObject("hDiffSideBand"); | |
1132 | TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb"); | |
1133 | ||
1134 | Double_t deltaM=candDstar->DeltaInvMass(); | |
1135 | //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar); | |
1136 | Double_t ptD=candDstar->Pt(); | |
1137 | Double_t ptjet=jet->Pt(); | |
1138 | Double_t dPhi=jet->Phi()-candDstar->Phi(); | |
1139 | Double_t deltaR=DeltaR(candDstar,jet); | |
1140 | if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi(); | |
1141 | if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi(); | |
1142 | ||
1143 | Int_t isSideBand=IsDzeroSideBand(candDstar); | |
1144 | //fill the background histogram with the side bands or when looking at MC with the real background | |
1145 | if(isSideBand==1){ | |
1146 | hDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background | |
1147 | //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi); | |
1148 | ||
1149 | if(deltaR<fJetRadius){ // evaluate in the near side | |
1150 | //hzptDB->Fill(Z(candDstar,jet),deltaM,ptD); | |
1151 | hPtJetWithDsb->Fill(ptjet,deltaM,ptD); | |
1152 | } | |
1153 | } | |
1154 | } | |
1155 | ||
1156 | //_______________________________________________________________________________ | |
1157 | ||
1158 | void AliAnalysisTaskFlavourFilterAndJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg, AliEmcalJet *jet){ | |
1159 | ||
1160 | //need mass, deltaR, pt jet, ptD | |
1161 | //two cases: D0 and Dstar | |
1162 | ||
1163 | Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent()); | |
1164 | if(!isselected) return; | |
1165 | ||
1166 | Double_t ptD=candbg->Pt(); | |
1167 | Double_t ptjet=jet->Pt(); | |
1168 | Double_t deltaR=DeltaR(candbg,jet); | |
1169 | ||
1170 | TH2F* hInvMassptDbg=(TH2F*)fmyOutput->FindObject("hInvMassptDbg"); | |
1171 | TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb"); | |
1172 | ||
1173 | ||
1174 | if(fCandidateType==kDstartoKpipi){ | |
1175 | AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg; | |
1176 | Double_t deltaM=dstarbg->DeltaInvMass(); | |
1177 | hInvMassptDbg->Fill(deltaM,ptD); | |
1178 | if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,deltaM,ptD); | |
1179 | } | |
1180 | ||
1181 | if(fCandidateType==kD0toKpi){ | |
1182 | Double_t masses[2]={0.,0.}; | |
1183 | Int_t pdgdaughtersD0[2]={211,321};//pi,K | |
1184 | Int_t pdgdaughtersD0bar[2]={321,211};//K,pi | |
1185 | ||
1186 | masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0 | |
1187 | masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar | |
1188 | ||
1189 | ||
1190 | if(isselected==1 || isselected==3) { | |
1191 | if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[0],ptD); | |
1192 | hInvMassptDbg->Fill(masses[0],ptD); | |
1193 | } | |
1194 | if(isselected>=2) { | |
1195 | if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[1],ptD); | |
1196 | hInvMassptDbg->Fill(masses[1],ptD); | |
1197 | } | |
1198 | ||
1199 | ||
1200 | } | |
1201 | } | |
1202 | ||
1203 | //_______________________________________________________________________________ | |
1204 | ||
1205 | Float_t AliAnalysisTaskFlavourFilterAndJetCorrelations::DeltaR(AliVParticle *p1, AliVParticle *p2) const { | |
1206 | //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2) | |
1207 | ||
1208 | if(!p1 || !p2) return -1; | |
1209 | Double_t phi1=p1->Phi(),eta1=p1->Eta(); | |
1210 | Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ; | |
1211 | ||
1212 | Double_t dPhi=phi1-phi2; | |
1213 | if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi()); | |
1214 | if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi()); | |
1215 | ||
1216 | Double_t dEta=eta1-eta2; | |
1217 | Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi ); | |
1218 | return deltaR; | |
1219 | ||
1220 | } | |
1221 | ||
1222 | //_______________________________________________________________________________ | |
1223 | ||
1224 | Int_t AliAnalysisTaskFlavourFilterAndJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){ | |
1225 | ||
1226 | Double_t ptD=candDstar->Pt(); | |
1227 | Int_t bin = fCuts->PtBin(ptD); | |
1228 | if (bin < 0){ | |
1229 | // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config. | |
1230 | bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this? | |
1231 | return -1; // out of bounds | |
1232 | } | |
1233 | ||
1234 | Double_t invM=candDstar->InvMassD0(); | |
1235 | Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass(); | |
1236 | ||
1237 | Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin]; | |
1238 | Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin]; | |
1239 | ||
1240 | if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1; | |
1241 | else return 0; | |
1242 | ||
1243 | } |