]>
Commit | Line | Data |
---|---|---|
954ac830 | 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 in Jets Analysis | |
18 | // | |
19 | // The D* (+ and -) is reconstructed inside jets. Two different cuts are | |
20 | // implemented: | |
21 | // | |
22 | // 1) C.Ivan D* cuts adapted for correlation analysis | |
23 | // 2) Topological cut enforcing the correlation D0-softPion pt + relaxed | |
24 | // CosThetaStar. This second should be better for correlations. | |
25 | // | |
26 | // USAGE: | |
27 | // | |
28 | // The analysis is performed separately for D*+ and D*-. A flag in the .C is | |
29 | // taking care to decide which analysis. | |
30 | // | |
31 | // The cut number 2 can be activaded with a flag in the .C (not active in this version) | |
32 | // Cuts 1) are the default. | |
33 | // | |
34 | // The task requires reconstructed jets in the AODs | |
35 | // | |
36 | //----------------------------------------------------------------------- | |
37 | // Author A.Grelli | |
38 | // ERC-QGP Utrecht University - a.grelli@uu.nl | |
39 | //----------------------------------------------------------------------- | |
40 | ||
41 | #include <TDatabasePDG.h> | |
42 | #include <TParticle.h> | |
43 | #include <TVector3.h> | |
44 | ||
45 | #include "AliAnalysisTaskSEDStarJets.h" | |
46 | #include "AliStack.h" | |
47 | #include "AliMCEvent.h" | |
48 | #include "AliAODMCHeader.h" | |
49 | #include "AliAnalysisManager.h" | |
50 | #include "AliAODHandler.h" | |
51 | #include "AliLog.h" | |
52 | #include "AliAODVertex.h" | |
53 | #include "AliAODJet.h" | |
54 | #include "AliAODRecoDecay.h" | |
55 | #include "AliAODRecoDecayHF.h" | |
56 | #include "AliAODRecoDecayHF2Prong.h" | |
57 | #include "AliESDtrack.h" | |
58 | #include "AliAODMCParticle.h" | |
59 | ||
60 | ClassImp(AliAnalysisTaskSEDStarJets) | |
61 | ||
62 | //__________________________________________________________________________ | |
63 | AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets() : | |
64 | AliAnalysisTaskSE(), | |
954ac830 | 65 | fCountReco(0), |
66 | fCountRecoAcc(0), | |
67 | fCountRecoITSClusters(0), | |
68 | fCountRecoPPR(0), | |
69 | fCountDStar(0), | |
954ac830 | 70 | fEvents(0), |
71 | fMinITSClusters(0), | |
72 | fComputeD0(kTRUE), | |
feaae220 | 73 | fUseMCInfo(kTRUE), |
954ac830 | 74 | ftopologicalCut(kFALSE), |
75 | fRequireNormalization(kTRUE), | |
76 | fLorentzTrack1(0,0,0,0), | |
77 | fLorentzTrack2(0,0,0,0), | |
78 | fLorentzTrack3(0,0,0,0), | |
79 | fLorentzTrack4(0,0,0,0), | |
80 | fOutput(0), | |
81 | fD0ptvsSoftPtSignal(0), | |
82 | fD0ptvsSoftPt(0), | |
83 | ftrigger(0), | |
84 | fPtPion(0), | |
85 | fInvMass(0), | |
86 | fRECOPtDStar(0), | |
87 | fDStar(0), | |
88 | fDiff(0), | |
89 | fDiffSideBand(0), | |
90 | fDStarMass(0), | |
91 | fPhi(0), | |
92 | fPhiBkg(0), | |
93 | fTrueDiff(0), | |
94 | fResZ(0), | |
95 | fResZBkg(0), | |
96 | fcharmpt(0), | |
97 | fdstarE(0), | |
98 | fEjet(0), | |
99 | fPhijet(0), | |
100 | fEtaJet(0), | |
c088a94c | 101 | fdstarpt(0) |
954ac830 | 102 | { |
103 | // | |
104 | // Default ctor | |
105 | // | |
106 | } | |
107 | //___________________________________________________________________________ | |
108 | AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets(const Char_t* name) : | |
109 | AliAnalysisTaskSE(name), | |
954ac830 | 110 | fCountReco(0), |
111 | fCountRecoAcc(0), | |
112 | fCountRecoITSClusters(0), | |
113 | fCountRecoPPR(0), | |
114 | fCountDStar(0), | |
954ac830 | 115 | fEvents(0), |
116 | fMinITSClusters(0), | |
117 | fComputeD0(kTRUE), | |
feaae220 | 118 | fUseMCInfo(kTRUE), |
954ac830 | 119 | ftopologicalCut(kFALSE), |
120 | fRequireNormalization(kTRUE), | |
121 | fLorentzTrack1(0,0,0,0), | |
122 | fLorentzTrack2(0,0,0,0), | |
123 | fLorentzTrack3(0,0,0,0), | |
124 | fLorentzTrack4(0,0,0,0), | |
125 | fOutput(0), | |
126 | fD0ptvsSoftPtSignal(0), | |
127 | fD0ptvsSoftPt(0), | |
128 | ftrigger(0), | |
129 | fPtPion(0), | |
130 | fInvMass(0), | |
131 | fRECOPtDStar(0), | |
132 | fDStar(0), | |
133 | fDiff(0), | |
134 | fDiffSideBand(0), | |
135 | fDStarMass(0), | |
136 | fPhi(0), | |
137 | fPhiBkg(0), | |
138 | fTrueDiff(0), | |
139 | fResZ(0), | |
140 | fResZBkg(0), | |
141 | fcharmpt(0), | |
142 | fdstarE(0), | |
143 | fEjet(0), | |
144 | fPhijet(0), | |
145 | fEtaJet(0), | |
c088a94c | 146 | fdstarpt(0) |
954ac830 | 147 | { |
148 | // | |
149 | // Constructor. Initialization of Inputs and Outputs | |
150 | // | |
151 | Info("AliAnalysisTaskSEDStarJets","Calling Constructor"); | |
152 | ||
153 | DefineOutput(1,TList::Class()); | |
154 | } | |
155 | ||
156 | //___________________________________________________________________________ | |
157 | AliAnalysisTaskSEDStarJets& AliAnalysisTaskSEDStarJets::operator=(const AliAnalysisTaskSEDStarJets& c) | |
158 | { | |
159 | // | |
160 | // Assignment operator | |
161 | // | |
162 | if (this!=&c) { | |
163 | AliAnalysisTaskSE::operator=(c) ; | |
164 | } | |
165 | ||
166 | return *this; | |
167 | } | |
168 | ||
169 | //___________________________________________________________________________ | |
170 | AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets(const AliAnalysisTaskSEDStarJets& c) : | |
171 | AliAnalysisTaskSE(c), | |
954ac830 | 172 | fCountReco(c.fCountReco), |
173 | fCountRecoAcc(c.fCountRecoAcc), | |
174 | fCountRecoITSClusters(c.fCountRecoITSClusters), | |
175 | fCountRecoPPR(c.fCountRecoPPR), | |
176 | fCountDStar(c.fCountDStar), | |
954ac830 | 177 | fEvents(c.fEvents), |
178 | fMinITSClusters(c.fMinITSClusters), | |
179 | fComputeD0(c.fComputeD0), | |
feaae220 | 180 | fUseMCInfo(c.fUseMCInfo), |
954ac830 | 181 | ftopologicalCut(c.ftopologicalCut), |
182 | fRequireNormalization(c.fRequireNormalization), | |
183 | fLorentzTrack1(c.fLorentzTrack1), | |
184 | fLorentzTrack2(c.fLorentzTrack2), | |
185 | fLorentzTrack3(c.fLorentzTrack3), | |
186 | fLorentzTrack4(c.fLorentzTrack4), | |
187 | fOutput(c.fOutput), | |
188 | fD0ptvsSoftPtSignal(c.fD0ptvsSoftPtSignal), | |
189 | fD0ptvsSoftPt(c.fD0ptvsSoftPt), | |
190 | ftrigger(c.ftrigger), | |
191 | fPtPion(c.fPtPion), | |
192 | fInvMass(c.fInvMass), | |
193 | fRECOPtDStar(c.fRECOPtDStar), | |
194 | fDStar(c.fDStar), | |
195 | fDiff(c.fDiff), | |
196 | fDiffSideBand(c.fDiffSideBand), | |
197 | fDStarMass(c.fDStarMass), | |
198 | fPhi(c.fPhi), | |
199 | fPhiBkg(c.fPhiBkg), | |
200 | fTrueDiff(c.fTrueDiff), | |
201 | fResZ(c.fResZ), | |
202 | fResZBkg(c.fResZBkg), | |
203 | fcharmpt(c.fcharmpt), | |
204 | fdstarE(c.fdstarE), | |
205 | fEjet(c.fEjet), | |
206 | fPhijet(c.fPhijet), | |
207 | fEtaJet(c.fEtaJet), | |
c088a94c | 208 | fdstarpt(c.fdstarpt) |
954ac830 | 209 | |
210 | { | |
211 | // | |
212 | // Copy Constructor | |
213 | // | |
214 | } | |
215 | ||
216 | //___________________________________________________________________________ | |
217 | AliAnalysisTaskSEDStarJets::~AliAnalysisTaskSEDStarJets() { | |
218 | // destructor | |
219 | ||
220 | Info("~AliAnalysisTaskSEDStarJets","Calling Destructor"); | |
221 | if (fOutput) { | |
222 | delete fOutput; | |
223 | fOutput = 0; | |
224 | } | |
225 | } | |
226 | ||
227 | //_________________________________________________ | |
228 | void AliAnalysisTaskSEDStarJets::UserExec(Option_t *) | |
229 | { | |
230 | // user exec | |
231 | if (!fInputEvent) { | |
232 | Error("UserExec","NO EVENT FOUND!"); | |
233 | return; | |
234 | } | |
235 | ||
236 | // Load the event | |
954ac830 | 237 | AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent); |
238 | ||
239 | TClonesArray *arrayVerticesHF=0; | |
240 | ||
241 | if(!aodEvent && AODEvent() && IsStandardAOD()) { | |
242 | // In case there is an AOD handler writing a standard AOD, use the AOD | |
243 | // event in memory rather than the input (ESD) event. | |
244 | aodEvent = dynamic_cast<AliAODEvent*> (AODEvent()); | |
245 | // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root) | |
246 | // have to taken from the AOD event hold by the AliAODExtension | |
247 | AliAODHandler* aodHandler = (AliAODHandler*) | |
248 | ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler()); | |
249 | if(aodHandler->GetExtensions()) { | |
250 | AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root"); | |
251 | AliAODEvent *aodFromExt = ext->GetAOD(); | |
252 | arrayVerticesHF=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi"); | |
253 | } | |
254 | } else { | |
255 | arrayVerticesHF=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi"); | |
256 | } | |
257 | ||
258 | if (!arrayVerticesHF){ | |
259 | AliInfo("Could not find array of HF vertices, skipping the event"); | |
260 | return; | |
261 | }else AliDebug(2, Form("Found %d vertices",arrayVerticesHF->GetEntriesFast())); | |
262 | ||
7c23877d | 263 | // fix for temporary bug in ESDfilter |
264 | // the AODs with null vertex pointer didn't pass the PhysSel | |
265 | if(!aodEvent->GetPrimaryVertex()) return; | |
266 | ||
267 | fEvents++; | |
268 | AliDebug(2,Form("Event %d",fEvents)); | |
269 | if (fEvents%10000 ==0) AliDebug(2,Form("Event %d",fEvents)); | |
270 | ||
271 | ||
954ac830 | 272 | // Simulate a jet triggered sample |
273 | TClonesArray *arrayofJets = (TClonesArray*)aodEvent->GetJets(); | |
274 | if(aodEvent->GetNJets()<=0) return; | |
275 | AliInfo("found a jet: processing D* in jet analysis"); | |
954ac830 | 276 | |
c088a94c | 277 | // AOD primary vertex |
278 | AliAODVertex *prVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex(); | |
279 | Double_t primaryPos[3]; | |
280 | prVtx->GetXYZ(primaryPos); | |
281 | Bool_t vtxFlag = kTRUE; | |
282 | TString title=prVtx->GetTitle(); | |
283 | if(!title.Contains("VertexerTracks")) vtxFlag=kFALSE; | |
284 | ||
954ac830 | 285 | // counters for efficiencies |
954ac830 | 286 | Int_t icountReco = 0; |
287 | Int_t icountRecoAcc = 0; | |
288 | Int_t icountRecoITSClusters = 0; | |
289 | Int_t icountRecoPPR = 0; | |
290 | Int_t fiDstar = 0; | |
291 | Int_t fDStarD0 = 0; | |
feaae220 | 292 | |
293 | // TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName())); | |
294 | ||
295 | if(fUseMCInfo){ | |
296 | TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName())); | |
297 | //loop on the MC event - some basic MC info on D*, D0 and soft pion | |
298 | if (!mcArray) AliError("Could not find Monte-Carlo in AOD"); | |
954ac830 | 299 | |
feaae220 | 300 | for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) { |
301 | AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart)); | |
302 | if (!mcPart) { | |
303 | AliWarning("Particle not found in tree, skipping"); | |
304 | continue; | |
305 | } | |
c088a94c | 306 | |
feaae220 | 307 | // charm pt |
308 | if(TMath::Abs(mcPart->GetPdgCode())==4){ | |
309 | fcharmpt->Fill(mcPart->Pt()); | |
c088a94c | 310 | } |
311 | ||
feaae220 | 312 | // fill energy and pt for D* in acceptance with correct prongs |
313 | Bool_t isOk = DstarInMC(mcPart,mcArray); | |
c088a94c | 314 | |
feaae220 | 315 | if (isOk){ //D* |
316 | AliDebug(2, "Found a DStar in MC with correct prongs and in acceptance"); | |
317 | fdstarE ->Fill(mcPart->E()); | |
318 | fdstarpt->Fill(mcPart->Pt()); | |
954ac830 | 319 | |
feaae220 | 320 | // check the MC-Acceptance level cuts |
321 | // since standard CF functions are not applicable, using Kine Cuts on daughters | |
322 | ||
323 | Int_t daughter0 = mcPart->GetDaughter(0); | |
324 | Int_t daughter1 = mcPart->GetDaughter(1); | |
325 | ||
326 | AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1)); | |
327 | ||
328 | AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0)); | |
329 | AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1)); | |
330 | ||
331 | Double_t eta0 = mcPartDaughter0->Eta(); | |
332 | Double_t eta1 = mcPartDaughter1->Eta(); | |
333 | Double_t y0 = mcPartDaughter0->Y(); | |
334 | Double_t y1 = mcPartDaughter1->Y(); | |
335 | Double_t pt0 = mcPartDaughter0->Pt(); | |
336 | Double_t pt1 = mcPartDaughter1->Pt(); | |
337 | ||
338 | AliDebug(2, Form("D* Daughter 0: eta = %f, y = %f, pt = %f", eta0, y0, pt0)); | |
339 | AliDebug(2, Form("D* Daughter 1: eta = %f, y = %f, pt = %f", eta1, y1, pt1)); | |
340 | ||
341 | Int_t daughD00 = 0; | |
342 | Int_t daughD01 = 0; | |
343 | ||
344 | // D0 daughters - do not need to check D0-kpi, already done | |
345 | ||
346 | if(TMath::Abs(mcPartDaughter0->GetPdgCode())==421){ | |
347 | daughD00 = mcPartDaughter0->GetDaughter(0); | |
348 | daughD01 = mcPartDaughter0->GetDaughter(1); | |
349 | }else{ | |
350 | daughD00 = mcPartDaughter1->GetDaughter(0); | |
351 | daughD01 = mcPartDaughter1->GetDaughter(1); | |
352 | } | |
353 | ||
354 | AliAODMCParticle* mcD0PartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD00)); | |
355 | AliAODMCParticle* mcD0PartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD01)); | |
356 | ||
357 | if (!mcD0PartDaughter0 || !mcD0PartDaughter1) { | |
358 | AliWarning("At least one Daughter Particle not found in tree, but it should be, this check was already done..."); | |
359 | } | |
360 | ||
361 | // D0 daughters - needed for acceptance | |
362 | Double_t pD0pt0 = mcD0PartDaughter0->Pt(); | |
363 | Double_t pD0pt1 = mcD0PartDaughter1->Pt(); | |
364 | Double_t pD0eta0 = mcD0PartDaughter0->Eta(); | |
365 | Double_t pD0eta1 = mcD0PartDaughter1->Eta(); | |
366 | ||
367 | // ACCEPTANCE REQUESTS --------- | |
368 | ||
369 | // soft pion | |
370 | Bool_t daught1inAcceptance = (TMath::Abs(eta1) <= 0.9 && pt1 > 0.05); | |
371 | // Do daughters | |
372 | Bool_t theD0daught0inAcceptance = (TMath::Abs(pD0eta0) <= 0.9 && pD0pt0 >= 0.1); | |
373 | Bool_t theD0daught1inAcceptance = (TMath::Abs(pD0eta1) <= 0.9 && pD0pt1 >= 0.1); | |
374 | ||
375 | if (daught1inAcceptance && theD0daught0inAcceptance && theD0daught1inAcceptance) { | |
376 | ||
377 | AliDebug(2, "Daughter particles in acceptance"); | |
378 | ||
379 | // check on the vertex | |
380 | if (vtxFlag){ | |
381 | printf("Vertex cut passed 2\n"); | |
382 | fDStarD0++; | |
383 | Bool_t refitFlag = kTRUE; | |
384 | for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){ | |
385 | AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod); | |
386 | ||
387 | // refit only for D0 daughters | |
388 | if ((track->GetLabel() == daughD00) || (track->GetLabel() == daughD01)) { | |
389 | if(!(track->GetStatus()&AliESDtrack::kTPCrefit) || !(track->GetStatus()&AliESDtrack::kITSrefit)) { | |
390 | refitFlag = kFALSE; | |
391 | } | |
c088a94c | 392 | } |
393 | } | |
feaae220 | 394 | if (refitFlag){ |
395 | printf("Refit cut passed\n"); | |
396 | } | |
397 | else{ | |
398 | AliDebug(3,"Refit cut not passed\n"); | |
399 | } | |
c088a94c | 400 | } |
401 | else{ | |
feaae220 | 402 | AliDebug(3,"Vertex cut not passed\n"); |
403 | } | |
954ac830 | 404 | } |
c088a94c | 405 | else{ |
feaae220 | 406 | AliDebug(3,"Acceptance cut not passed\n"); |
407 | } | |
954ac830 | 408 | } |
954ac830 | 409 | } |
feaae220 | 410 | |
411 | AliDebug(2, Form("Found %i MC particles that are D* in Kpipi and satisfy Acc cuts!!",fDStarD0)); | |
412 | fCountDStar += fDStarD0; | |
413 | ||
414 | } // End of MC | |
954ac830 | 415 | |
416 | // Now perform the D* in jet reconstruction | |
feaae220 | 417 | |
954ac830 | 418 | // Normalization factor |
419 | if(fRequireNormalization){ | |
420 | ftrigger->Fill(1); | |
421 | } | |
422 | ||
423 | Int_t nJets = 0; // for reco D0 countings | |
424 | ||
425 | for (Int_t iJets = 0; iJets<arrayofJets->GetEntriesFast(); iJets++) { // main loop on jets | |
426 | ||
427 | AliAODJet* jet = (AliAODJet*)arrayofJets->At(iJets); | |
428 | ||
429 | //jets variables | |
430 | Double_t ejet = jet->E(); | |
431 | Double_t phiJet = jet->Phi(); | |
432 | Double_t etaJet = jet->Eta(); | |
433 | ||
434 | // fill energy, eta and phi of the jet | |
435 | fEjet ->Fill(ejet); | |
436 | fPhijet ->Fill(phiJet); | |
437 | fEtaJet ->Fill(etaJet); | |
438 | ||
439 | nJets++; | |
440 | ||
441 | // loop over the tracks to search for candidates soft pions | |
442 | for (Int_t i=0; i<aodEvent->GetNTracks(); i++){ | |
443 | ||
444 | AliAODTrack* aodTrack = aodEvent->GetTrack(i); | |
954ac830 | 445 | Double_t vD0[4] = {0.,0.,0.,0.}; |
446 | Double_t vD0bar[4] ={0.,0.,0.,0.}; | |
447 | ||
448 | Int_t pCharge = aodTrack->Charge(); | |
449 | ||
450 | // few selections on soft pion | |
451 | Bool_t tPrimCand = aodTrack->IsPrimaryCandidate(); // is it primary? | |
452 | ||
453 | if(aodTrack->Pt()>= 5 || aodTrack->Pt()< 0.08 ) continue; //cut on soft pion pt VERY large | |
454 | //~ D*s of pt >80GeV with a soft pion of 5GeV! | |
455 | if(TMath::Abs(aodTrack->Eta())>0.9) continue; | |
456 | ||
c088a94c | 457 | // if D*+ analysis tha D0 and pi+ otherwise pi- |
954ac830 | 458 | if(fComputeD0 && pCharge!= 1 ) continue; |
954ac830 | 459 | if(!fComputeD0 && pCharge!= -1 ) continue; |
460 | ||
461 | if(tPrimCand && arrayVerticesHF->GetEntriesFast()>0){ // isPion and is Primary, no PID for now | |
462 | ||
463 | // label to the candidate soft pion | |
feaae220 | 464 | Int_t pLabel = -1; |
465 | if(fUseMCInfo) pLabel = aodTrack->GetLabel(); | |
954ac830 | 466 | |
467 | // prepare the TLorentz vector for the pion | |
468 | Float_t pionMass = TDatabasePDG::Instance()->GetParticle(211)->Mass(); | |
469 | fLorentzTrack3.SetPxPyPzE(aodTrack->Px(),aodTrack->Py(),aodTrack->Pz(),aodTrack->E(pionMass)); | |
470 | ||
471 | // search for the D0 | |
472 | for (Int_t iVertex = 0; iVertex<arrayVerticesHF->GetEntriesFast(); iVertex++) { | |
473 | ||
474 | Double_t invM = 0; | |
475 | Double_t invMDStar = 0; | |
c088a94c | 476 | Double_t dPhi = 0; |
954ac830 | 477 | |
478 | AliAODRecoDecayHF2Prong* vtx = (AliAODRecoDecayHF2Prong*)arrayVerticesHF->At(iVertex); | |
479 | ||
480 | Double_t pt = vtx->Pt(); | |
c088a94c | 481 | |
482 | // acceptance variables | |
483 | Bool_t acceptanceProng0 = (TMath::Abs(vtx->EtaProng(0))<= 0.9 && vtx->PtProng(0) >= 0.1); | |
484 | Bool_t acceptanceProng1 = (TMath::Abs(vtx->EtaProng(1))<= 0.9 && vtx->PtProng(1) >= 0.1); | |
485 | Bool_t acceptanceSoftPi = (TMath::Abs(aodTrack->Eta())<= 0.9 && aodTrack->Pt() >= 0.05); | |
954ac830 | 486 | |
954ac830 | 487 | Int_t pdgDgD0toKpi[2]={321,211}; |
c088a94c | 488 | |
489 | Int_t mcLabel =-1; | |
c088a94c | 490 | |
491 | Bool_t isDStar = kFALSE; // just to count | |
492 | ||
feaae220 | 493 | //switch of MC for real data |
494 | if(fUseMCInfo){ | |
495 | TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName())); | |
496 | mcLabel = vtx->MatchToMC(421, mcArray,2,pdgDgD0toKpi) ; //MC D0 | |
497 | // matching to MC D* | |
498 | if(mcLabel !=-1 && pLabel!=-1 && nJets ==1) { // count only once in case of multijets | |
499 | ||
500 | // search for a D0 and a pi with mother and check it is a D* | |
501 | AliAODMCParticle* theMCpion = (AliAODMCParticle*)mcArray->At(pLabel); | |
502 | Int_t motherMCPion = theMCpion->GetMother(); | |
503 | AliAODMCParticle* theMCD0 = (AliAODMCParticle*)mcArray->At(mcLabel); | |
504 | Int_t motherMCD0 = theMCD0->GetMother(); | |
505 | ||
506 | if(motherMCPion!=-1 && (motherMCD0 == motherMCPion)){ | |
507 | AliAODMCParticle* mcMother = (AliAODMCParticle*)mcArray->At(motherMCPion); | |
508 | if(TMath::Abs(mcMother->GetPdgCode()) == 413 && TMath::Abs(theMCpion->GetPdgCode())==211){ | |
509 | isDStar = kTRUE; | |
510 | fiDstar++; | |
511 | } | |
c088a94c | 512 | } |
513 | } | |
514 | } | |
954ac830 | 515 | |
feaae220 | 516 | |
c088a94c | 517 | if (acceptanceProng0 && acceptanceProng1 && acceptanceSoftPi) { |
954ac830 | 518 | |
519 | AliDebug(2,"D0 reco daughters in acceptance"); | |
c088a94c | 520 | if(isDStar && nJets==1) icountRecoAcc++; |
954ac830 | 521 | |
c088a94c | 522 | // D0 daughters |
954ac830 | 523 | AliAODTrack *track0 = (AliAODTrack*)vtx->GetDaughter(0); |
524 | AliAODTrack *track1 = (AliAODTrack*)vtx->GetDaughter(1); | |
525 | ||
526 | // check for ITS refit (already required at the ESD filter level ) | |
c088a94c | 527 | Bool_t kRefitITS = kTRUE; |
954ac830 | 528 | |
529 | if((!(track0->GetStatus()&AliESDtrack::kITSrefit)|| (!(track1->GetStatus()&AliESDtrack::kITSrefit)))) { | |
530 | kRefitITS = kFALSE; | |
c088a94c | 531 | } |
954ac830 | 532 | |
c088a94c | 533 | Int_t ncls0=0,ncls1=0,ncls2=0; |
534 | ||
954ac830 | 535 | for(Int_t l=0;l<6;l++) { |
536 | if(TESTBIT(track0->GetITSClusterMap(),l)) ncls0++; | |
537 | if(TESTBIT(track1->GetITSClusterMap(),l)) ncls1++; | |
c088a94c | 538 | if(TESTBIT(aodTrack->GetITSClusterMap(),l)) ncls2++; |
539 | } | |
540 | ||
541 | // clusters in ITS for D0 daugthers and soft pion | |
feaae220 | 542 | if (ncls0 >= fMinITSClusters && ncls1 >= fMinITSClusters && ncls2>=4 && kRefitITS) { |
954ac830 | 543 | |
c088a94c | 544 | if(isDStar && nJets==1) icountRecoITSClusters++; |
954ac830 | 545 | |
546 | // D0 cutting varibles | |
547 | Double_t cosThetaStar = 9999.; | |
548 | Double_t pTpi = 0.; | |
549 | Double_t pTK = 0.; | |
c088a94c | 550 | Double_t d01 = 0; |
551 | Double_t d00 = 0; | |
954ac830 | 552 | |
553 | // D0, D0bar | |
c088a94c | 554 | if (fComputeD0){ |
555 | cosThetaStar = vtx->CosThetaStarD0(); | |
954ac830 | 556 | pTpi = vtx->PtProng(0); |
c088a94c | 557 | pTK = vtx->PtProng(1); |
558 | d01 = vtx->Getd0Prong(0)*1E4; | |
559 | d00 = vtx->Getd0Prong(1)*1E4; | |
560 | }else{ | |
561 | cosThetaStar = vtx->CosThetaStarD0bar(); | |
954ac830 | 562 | pTpi = vtx->PtProng(1); |
c088a94c | 563 | pTK = vtx->PtProng(0); |
564 | d01 = vtx->Getd0Prong(1)*1E4; | |
565 | d00 = vtx->Getd0Prong(0)*1E4; | |
954ac830 | 566 | } |
c088a94c | 567 | |
954ac830 | 568 | AliDebug(2,"D0 reco daughters in acceptance"); |
c088a94c | 569 | |
954ac830 | 570 | Double_t dca = vtx->GetDCA()*1E4; |
954ac830 | 571 | Double_t d0d0 = vtx->Prodd0d0()*1E8; |
c088a94c | 572 | |
954ac830 | 573 | TVector3 mom(vtx->Px(),vtx->Py(),vtx->Pz()); |
574 | TVector3 flight((vtx->Xv())- primaryPos[0],(vtx->Yv())- primaryPos[1],(vtx->Zv())- primaryPos[2]); | |
575 | Double_t pta = mom.Angle(flight); | |
576 | Double_t cosPointingAngle = TMath::Cos(pta); | |
c088a94c | 577 | |
578 | // Crstian Ivan D* cuts. Multidimensional optimization | |
579 | Double_t cuts[7] = {9999999., 1.1, 0., 9999999.,999999., 9999999., 0.}; | |
954ac830 | 580 | |
581 | if (pt <= 1){ // first bin not optimized | |
582 | cuts[0] = 400; | |
583 | cuts[1] = 0.8; | |
584 | cuts[2] = 0.21; | |
585 | cuts[3] = 500; | |
c088a94c | 586 | cuts[4] = 500; |
587 | cuts[5] = -20000; | |
588 | cuts[6] = 0.6; | |
954ac830 | 589 | } |
590 | else if (pt > 1 && pt <= 2){ | |
c088a94c | 591 | cuts[0] = 200; |
954ac830 | 592 | cuts[1] = 0.7; |
593 | cuts[2] = 0.8; | |
c088a94c | 594 | cuts[3] = 210; |
595 | cuts[4] = 210; | |
596 | cuts[5] = -20000; | |
597 | cuts[6] = 0.9; | |
954ac830 | 598 | } |
599 | else if (pt > 2 && pt <= 3){ | |
600 | cuts[0] = 400; | |
601 | cuts[1] = 0.8; | |
602 | cuts[2] = 0.8; | |
c088a94c | 603 | cuts[3] = 420; |
604 | cuts[4] = 350; | |
605 | cuts[5] = -8500; | |
606 | cuts[6] = 0.9; | |
954ac830 | 607 | } |
608 | else if (pt > 3 && pt <= 5){ | |
c088a94c | 609 | cuts[0] = 160; |
610 | cuts[1] = 1.0; | |
954ac830 | 611 | cuts[2] = 1.2; |
c088a94c | 612 | cuts[3] = 420; |
613 | cuts[4] = 560; | |
614 | cuts[5] = -8500; | |
615 | cuts[6] = 0.9; | |
954ac830 | 616 | } |
617 | else if (pt > 5){ | |
618 | cuts[0] = 800; | |
619 | cuts[1] = 1.0; | |
620 | cuts[2] = 1.2; | |
c088a94c | 621 | cuts[3] = 700; |
622 | cuts[4] = 700; | |
623 | cuts[5] = 10000; | |
624 | cuts[6] = 0.9; | |
954ac830 | 625 | } |
626 | // apply D0 cuts | |
627 | ||
628 | if (dca < cuts[0] | |
c088a94c | 629 | && TMath::Abs(cosThetaStar) < cuts[1] |
630 | && pTpi > cuts[2] | |
631 | && pTK > cuts[2] | |
632 | && TMath::Abs(d01) < cuts[3] | |
633 | && TMath::Abs(d00) < cuts[4] | |
634 | && d0d0 < cuts[5] | |
635 | && cosPointingAngle > cuts[6] | |
636 | ){ | |
954ac830 | 637 | |
638 | if(fComputeD0){ // D0 from default | |
639 | ||
640 | if(vtx->ChargeProng(0)==-1){ | |
641 | fLorentzTrack1.SetPxPyPzE( vtx->PxProng(0),vtx->PyProng(0), vtx->PzProng(0),vtx->EProng(0,321) ); | |
642 | fLorentzTrack2.SetPxPyPzE( vtx->PxProng(1),vtx->PyProng(1), vtx->PzProng(1),vtx->EProng(1,211) ); | |
643 | }else{ | |
644 | fLorentzTrack1.SetPxPyPzE( vtx->PxProng(0),vtx->PyProng(0), vtx->PzProng(0),vtx->EProng(0,211) ); | |
645 | fLorentzTrack2.SetPxPyPzE( vtx->PxProng(1),vtx->PyProng(1), vtx->PzProng(1),vtx->EProng(1,321) ); | |
646 | } | |
647 | ||
954ac830 | 648 | vD0[0] = (fLorentzTrack1+fLorentzTrack2).Px(); |
649 | vD0[1] = (fLorentzTrack1+fLorentzTrack2).Py(); | |
650 | vD0[2] = (fLorentzTrack1+fLorentzTrack2).Pz(); | |
651 | vD0[3] = (fLorentzTrack1+fLorentzTrack2).E(); | |
652 | ||
653 | fLorentzTrack4.SetPxPyPzE(vD0[0],vD0[1],vD0[2],vD0[3]); | |
654 | ||
655 | }else{ // D0bar analysis | |
656 | ||
657 | if(vtx->ChargeProng(0)==-1){ | |
658 | fLorentzTrack1.SetPxPyPzE( vtx->PxProng(0),vtx->PyProng(0), vtx->PzProng(0),vtx->EProng(0,211) ); | |
659 | fLorentzTrack2.SetPxPyPzE( vtx->PxProng(1),vtx->PyProng(1), vtx->PzProng(1),vtx->EProng(1,321) ); | |
660 | }else{ | |
661 | fLorentzTrack1.SetPxPyPzE( vtx->PxProng(0),vtx->PyProng(0), vtx->PzProng(0),vtx->EProng(0,321) ); | |
662 | fLorentzTrack2.SetPxPyPzE( vtx->PxProng(1),vtx->PyProng(1), vtx->PzProng(1),vtx->EProng(1,211) ); | |
663 | } | |
664 | ||
954ac830 | 665 | vD0bar[0] = (fLorentzTrack1+fLorentzTrack2).Px(); |
666 | vD0bar[1] = (fLorentzTrack1+fLorentzTrack2).Py(); | |
667 | vD0bar[2] = (fLorentzTrack1+fLorentzTrack2).Pz(); | |
668 | vD0bar[3] = (fLorentzTrack1+fLorentzTrack2).E(); | |
669 | ||
670 | fLorentzTrack4.SetPxPyPzE(vD0bar[0],vD0bar[1],vD0bar[2],vD0bar[3]); | |
671 | } | |
672 | ||
673 | // D0-D0bar related variables | |
674 | ||
675 | invM = GetInvariantMass(fLorentzTrack1,fLorentzTrack2); | |
c088a94c | 676 | if(nJets==1) fInvMass->Fill(invM); |
954ac830 | 677 | |
c088a94c | 678 | if(isDStar && nJets==1) icountRecoPPR++; |
679 | ||
954ac830 | 680 | //DStar invariant mass |
681 | invMDStar = GetInvariantMassDStar(fLorentzTrack3,fLorentzTrack4); | |
682 | ||
683 | //conversion from phi TLorentzVerctor to phi aliroot | |
684 | Double_t kconvert = 0; | |
685 | Double_t phiDStar = (fLorentzTrack3 + fLorentzTrack4).Phi(); | |
686 | kconvert = phiDStar; | |
687 | if(phiDStar<0) kconvert = phiDStar + 2*(TMath::Pi()); | |
688 | phiDStar = kconvert; | |
689 | ||
690 | // dphi between jet and D* | |
691 | dPhi = phiJet - phiDStar; | |
692 | ||
693 | //Just for plotting pourposal | |
694 | if(dPhi<=-1.58) dPhi = TMath::Abs(dPhi); | |
695 | if(dPhi>4.72){ | |
696 | dPhi = dPhi-2*(TMath::Pi()); | |
697 | } | |
698 | ||
699 | //Alternative cutting procedure | |
700 | //the cut on cosThetaStar is to reduce near collinear | |
701 | //background from jet fragmentation | |
702 | ||
703 | Bool_t tCut = EvaluateCutOnPiD0pt(vtx,aodTrack); | |
704 | ||
705 | if(ftopologicalCut && tCut) continue; | |
706 | if(ftopologicalCut && TMath::Abs(cosThetaStar)>cuts[1]) continue; | |
707 | ||
708 | if(invM >= 1.829 && invM <= 1.901){ // ~4 sigma cut on D0 mass | |
709 | ||
c088a94c | 710 | if(isDStar && nJets==1) { |
711 | fDStarMass->Fill(invMDStar); | |
712 | fTrueDiff->Fill(invMDStar-invM); | |
954ac830 | 713 | } |
714 | ||
c088a94c | 715 | if(nJets==1) { // avoid double counting |
716 | fDStar->Fill(invMDStar); | |
717 | fDiff->Fill(invMDStar-invM); // M(Kpipi)-M(Kpi) | |
718 | } | |
719 | ||
954ac830 | 720 | // now the dphi signal and the fragmentation function |
feaae220 | 721 | if((invMDStar-invM)<=0.148 && (invMDStar-invM)>=0.142) { |
954ac830 | 722 | |
723 | //fill candidates D* and soft pion reco pt | |
c088a94c | 724 | if(nJets==1) fPtPion->Fill(aodTrack->Pt()); |
725 | if(nJets==1) fRECOPtDStar->Fill((fLorentzTrack3 + fLorentzTrack4).Pt()); | |
954ac830 | 726 | fPhi ->Fill(dPhi); |
727 | ||
728 | if(dPhi>=-0.5 && dPhi<=0.5){ // evaluate in the near side | |
729 | Double_t dStarMom = (fLorentzTrack3 + fLorentzTrack4).P(); | |
730 | Double_t zFrag = (TMath::Cos(dPhi)*dStarMom)/ejet; | |
731 | fResZ->Fill(TMath::Abs(zFrag)); | |
732 | } | |
733 | } | |
734 | } | |
735 | ||
736 | // evaluate side band background | |
feaae220 | 737 | SideBandBackground(invM, invMDStar, ejet, dPhi, nJets); |
954ac830 | 738 | |
739 | invM = 0; | |
740 | invMDStar = 0; | |
741 | ||
c088a94c | 742 | } // end PPR cuts |
743 | } // end ITS cluster | |
744 | } // end acceptance | |
745 | } // D0 cand | |
746 | } // softpion | |
747 | } // tracks | |
954ac830 | 748 | } // jets |
749 | ||
750 | AliDebug(2, Form("Found %i Reco particles that are D0!!",icountReco)); | |
751 | ||
752 | fCountReco+= fiDstar; | |
753 | fCountRecoAcc+= icountRecoAcc; | |
754 | fCountRecoITSClusters+= icountRecoITSClusters; | |
755 | fCountRecoPPR+= icountRecoPPR; | |
756 | ||
757 | PostData(1,fOutput); | |
758 | } | |
759 | ||
760 | //________________________________________ terminate ___________________________ | |
761 | void AliAnalysisTaskSEDStarJets::Terminate(Option_t*) | |
762 | { | |
763 | // The Terminate() function is the last function to be called during | |
764 | // a query. It always runs on the client, it can be used to present | |
765 | // the results graphically or save the results to file. | |
766 | ||
954ac830 | 767 | AliAnalysisTaskSE::Terminate(); |
768 | ||
954ac830 | 769 | AliInfo(Form("Found %i of that MC D*->D0pi(D0->kpi) are in acceptance, in %d events",fCountDStar,fEvents)); |
770 | ||
771 | AliInfo(Form("Found %i RECO particles after cuts that are DStar, in %d events",fCountReco,fEvents)); | |
c088a94c | 772 | AliInfo(Form("Among the above, found %i reco D* that are decaying in K+pi+pi and are in the requested acceptance, in %d events",fCountRecoAcc,fEvents)); |
773 | AliInfo(Form("Among the above, found %i reco D* that are decaying in K+pi+pi and have at least %d clusters in ITS, in %d events",fCountRecoITSClusters,fMinITSClusters,fEvents)); | |
774 | AliInfo(Form("Among the above, found %i reco D* that are decaying in K+pi+pi and satisfy D* cuts, in %d events",fCountRecoPPR,fEvents)); | |
954ac830 | 775 | |
776 | fOutput = dynamic_cast<TList*> (GetOutputData(1)); | |
777 | if (!fOutput) { | |
778 | printf("ERROR: fOutput not available\n"); | |
779 | return; | |
780 | } | |
781 | ||
954ac830 | 782 | fcharmpt = dynamic_cast<TH1F*>(fOutput->FindObject("fcharmpt")); |
783 | fdstarE = dynamic_cast<TH1F*>(fOutput->FindObject("fdstarE")); | |
784 | fdstarpt = dynamic_cast<TH1F*>(fOutput->FindObject("fdstarpt")); | |
785 | fDStarMass = dynamic_cast<TH1F*>(fOutput->FindObject("fDStarMass")); | |
786 | fTrueDiff = dynamic_cast<TH1F*>(fOutput->FindObject("fTrueDiff")); | |
787 | fInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass")); | |
788 | fPtPion = dynamic_cast<TH1F*>(fOutput->FindObject("fPtPion ")); | |
789 | fDStar = dynamic_cast<TH1F*>(fOutput->FindObject("fDStar")); | |
790 | fDiff = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff")); | |
791 | fDiffSideBand = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand")); | |
792 | ftrigger = dynamic_cast<TH1F*>(fOutput->FindObject("ftrigger")); | |
793 | fRECOPtDStar = dynamic_cast<TH1F*>(fOutput->FindObject("fRECOPtDStar")); | |
794 | fEjet = dynamic_cast<TH1F*>(fOutput->FindObject("fEjet")); | |
795 | fPhijet = dynamic_cast<TH1F*>(fOutput->FindObject("fPhijet")); | |
796 | fEtaJet = dynamic_cast<TH1F*>(fOutput->FindObject("fEtaJet")); | |
797 | fPhi = dynamic_cast<TH1F*>(fOutput->FindObject("fPhi")); | |
798 | fResZ = dynamic_cast<TH1F*>(fOutput->FindObject("fResZ")); | |
799 | fResZBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fResZBkg")); | |
800 | fPhiBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fPhiBkg")); | |
801 | fD0ptvsSoftPtSignal = dynamic_cast<TH2F*>(fOutput->FindObject("fD0ptvsSoftPtSignal")); | |
802 | fD0ptvsSoftPt = dynamic_cast<TH2F*>(fOutput->FindObject("fD0ptvsSoftPt")); | |
803 | ||
804 | } | |
805 | //___________________________________________________________________________ | |
806 | ||
807 | void AliAnalysisTaskSEDStarJets::UserCreateOutputObjects() { | |
808 | // output | |
809 | ||
810 | Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName()); | |
811 | ||
812 | //slot #1 | |
813 | OpenFile(1); | |
814 | fOutput = new TList(); | |
815 | fOutput->SetOwner(); | |
816 | // define histograms | |
817 | DefineHistoFroAnalysis(); | |
818 | ||
819 | return; | |
820 | } | |
954ac830 | 821 | |
c088a94c | 822 | //_______________________________D0 invariant mass________________________________ |
954ac830 | 823 | |
c088a94c | 824 | Double_t AliAnalysisTaskSEDStarJets::GetInvariantMass(TLorentzVector LorentzTrack1, TLorentzVector LorentzTrack2){ |
825 | ||
826 | return TMath::Abs((LorentzTrack1+LorentzTrack2).M()); // invariant mass of two tracks | |
827 | } | |
828 | //______________________________D* invariant mass_________________________________ | |
954ac830 | 829 | |
c088a94c | 830 | Double_t AliAnalysisTaskSEDStarJets::GetInvariantMassDStar(TLorentzVector LorentzTrack4, TLorentzVector LorentzTrack3){ |
831 | ||
832 | return TMath::Abs((LorentzTrack4+LorentzTrack3).M()); // invariant mass of two tracks | |
833 | } | |
834 | //_________________________________D* in MC _______________________________________ | |
954ac830 | 835 | |
c088a94c | 836 | Bool_t AliAnalysisTaskSEDStarJets::DstarInMC(AliAODMCParticle* const mcPart, TClonesArray* mcArray){ |
837 | // D* in MC | |
838 | ||
839 | Bool_t isOk = kFALSE; | |
840 | ||
841 | if(TMath::Abs(mcPart->GetPdgCode())!=413) return isOk; | |
954ac830 | 842 | |
843 | // getting the daughters | |
844 | Int_t daughter0 = mcPart->GetDaughter(0); | |
845 | Int_t daughter1 = mcPart->GetDaughter(1); | |
846 | ||
847 | AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1)); | |
848 | if (daughter0 == 0 || daughter1 == 0) { | |
c088a94c | 849 | AliDebug(2, "Error! the D* MC doesn't have correct daughters!!"); |
954ac830 | 850 | return isOk; |
851 | } | |
c088a94c | 852 | |
853 | if (TMath::Abs(daughter1 - daughter0) != 1) { // should be everytime true - see PDGdatabooklet | |
854 | AliDebug(2, "The D* MC doesn't come from a 2-prong decay, skipping!!"); | |
954ac830 | 855 | return isOk; |
856 | } | |
954ac830 | 857 | |
c088a94c | 858 | AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0)); |
859 | AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1)); | |
954ac830 | 860 | if (!mcPartDaughter0 || !mcPartDaughter1) { |
861 | AliWarning("At least one Daughter Particle not found in tree, skipping"); | |
862 | return isOk; | |
863 | } | |
864 | ||
c088a94c | 865 | if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==421 && |
866 | TMath::Abs(mcPartDaughter1->GetPdgCode())==211) && | |
867 | !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 && | |
868 | TMath::Abs(mcPartDaughter1->GetPdgCode())==421)) { | |
869 | AliDebug(2, "The D* MC doesn't come from a Kpi decay, skipping!!"); | |
870 | return isOk; | |
871 | } | |
872 | ||
873 | Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0 | |
874 | Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1 | |
875 | ||
876 | // getting vertex from daughters | |
877 | mcPartDaughter0->XvYvZv(vtx2daughter0); // cm | |
878 | mcPartDaughter1->XvYvZv(vtx2daughter1); //cm | |
879 | ||
880 | // check if the secondary vertex is the same for both | |
881 | if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) { | |
882 | AliError("Daughters have different secondary vertex, skipping the track"); | |
883 | return isOk; | |
884 | } | |
885 | ||
886 | AliAODMCParticle* neutralDaugh = mcPartDaughter0; | |
887 | ||
888 | if (!EvaluateIfD0toKpi(neutralDaugh,mcArray)) { | |
889 | AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!"); | |
890 | return isOk; | |
891 | } | |
954ac830 | 892 | |
893 | isOk = kTRUE; | |
894 | ||
895 | return isOk; | |
954ac830 | 896 | |
954ac830 | 897 | } |
954ac830 | 898 | |
c088a94c | 899 | Bool_t AliAnalysisTaskSEDStarJets::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, TClonesArray* mcArray)const{ |
954ac830 | 900 | |
c088a94c | 901 | // |
902 | // chack wether D0 is decaing into kpi | |
903 | // | |
954ac830 | 904 | |
c088a94c | 905 | Bool_t isHadronic = kFALSE; |
954ac830 | 906 | |
c088a94c | 907 | Int_t daughterD00 = neutralDaugh->GetDaughter(0); |
908 | Int_t daughterD01 = neutralDaugh->GetDaughter(1); | |
954ac830 | 909 | |
c088a94c | 910 | AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughterD00,daughterD01)); |
911 | if (daughterD00 == 0 || daughterD01 == 0) { | |
912 | AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!"); | |
913 | return isHadronic; | |
954ac830 | 914 | } |
c088a94c | 915 | |
916 | if (TMath::Abs(daughterD01 - daughterD00) != 1) { // should be everytime true - see PDGdatabooklet | |
917 | AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!"); | |
918 | return isHadronic; | |
954ac830 | 919 | } |
c088a94c | 920 | |
921 | AliAODMCParticle* mcPartDaughterD00 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD00)); | |
922 | AliAODMCParticle* mcPartDaughterD01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD01)); | |
923 | if (!mcPartDaughterD00 || !mcPartDaughterD01) { | |
954ac830 | 924 | AliWarning("At least one Daughter Particle not found in tree, skipping"); |
c088a94c | 925 | return isHadronic; |
926 | } | |
954ac830 | 927 | |
c088a94c | 928 | if (!(TMath::Abs(mcPartDaughterD00->GetPdgCode())==321 && |
929 | TMath::Abs(mcPartDaughterD01->GetPdgCode())==211) && | |
930 | !(TMath::Abs(mcPartDaughterD00->GetPdgCode())==211 && | |
931 | TMath::Abs(mcPartDaughterD01->GetPdgCode())==321)) { | |
932 | AliDebug(2, "The D* MC doesn't come from a Kpi decay, skipping!!"); | |
933 | return isHadronic; | |
934 | } | |
954ac830 | 935 | |
c088a94c | 936 | isHadronic = kTRUE; |
937 | ||
938 | ||
939 | return isHadronic; | |
954ac830 | 940 | |
941 | } | |
942 | ||
c088a94c | 943 | |
954ac830 | 944 | //___________________________________ hiostograms _______________________________________ |
945 | ||
946 | Bool_t AliAnalysisTaskSEDStarJets::DefineHistoFroAnalysis(){ | |
947 | ||
948 | // Invariant mass related histograms | |
949 | fInvMass = new TH1F("invMass","Kpi invariant mass distribution",1500,.5,3.5); | |
950 | fInvMass->SetStats(kTRUE); | |
951 | fInvMass->GetXaxis()->SetTitle("GeV/c"); | |
952 | fInvMass->GetYaxis()->SetTitle("Entries"); | |
953 | fOutput->Add(fInvMass); | |
954 | ||
955 | fDStar = new TH1F("invMassDStar","DStar invariant mass after D0 cuts ",600,1.8,2.4); | |
956 | fDStar->SetStats(kTRUE); | |
957 | fDStar->GetXaxis()->SetTitle("GeV/c"); | |
958 | fDStar->GetYaxis()->SetTitle("Entries"); | |
959 | fOutput->Add(fDStar); | |
960 | ||
961 | fDiff = new TH1F("Diff","M(kpipi)-M(kpi)",750,0.1,0.2); | |
962 | fDiff->SetStats(kTRUE); | |
963 | fDiff->GetXaxis()->SetTitle("M(kpipi)-M(kpi) GeV/c^2"); | |
964 | fDiff->GetYaxis()->SetTitle("Entries"); | |
965 | fOutput->Add(fDiff); | |
966 | ||
967 | fDiffSideBand = new TH1F("DiffSide","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2); | |
968 | fDiffSideBand->SetStats(kTRUE); | |
969 | fDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2"); | |
970 | fDiffSideBand->GetYaxis()->SetTitle("Entries"); | |
971 | fOutput->Add(fDiffSideBand); | |
972 | ||
973 | fDStarMass = new TH1F("RECODStar2","RECO DStar invariant mass distribution",750,1.5,2.5); | |
974 | fOutput->Add(fDStarMass); | |
975 | ||
976 | fTrueDiff = new TH1F("dstar","True Reco diff",750,0,0.2); | |
977 | fOutput->Add(fTrueDiff); | |
978 | ||
979 | // trigger normalization | |
980 | ftrigger = new TH1F("Normalization","Normalization factor for correlations",1,0,10); | |
981 | ftrigger->SetStats(kTRUE); | |
982 | fOutput->Add(ftrigger); | |
983 | ||
984 | //correlation fistograms | |
985 | fPhi = new TH1F("phi","Delta phi between Jet axis and DStar ",25,-1.57,4.72); | |
986 | fPhi->SetStats(kTRUE); | |
987 | fPhi->GetXaxis()->SetTitle("#Delta #phi (rad)"); | |
988 | fPhi->GetYaxis()->SetTitle("Entries"); | |
989 | fOutput->Add(fPhi); | |
990 | ||
991 | fPhiBkg = new TH1F("phiBkg","Delta phi between Jet axis and DStar background ",25,-1.57,4.72); | |
992 | fPhiBkg->SetStats(kTRUE); | |
993 | fPhiBkg->GetXaxis()->SetTitle("#Delta #phi (rad)"); | |
994 | fPhiBkg->GetYaxis()->SetTitle("Entries"); | |
995 | fOutput->Add(fPhiBkg); | |
996 | ||
997 | fRECOPtDStar = new TH1F("RECODStar1","RECO DStar pt distribution",600,0,15); | |
998 | fRECOPtDStar->SetStats(kTRUE); | |
999 | fRECOPtDStar->SetLineColor(2); | |
1000 | fRECOPtDStar->GetXaxis()->SetTitle("GeV/c"); | |
1001 | fRECOPtDStar->GetYaxis()->SetTitle("Entries"); | |
1002 | fOutput->Add(fRECOPtDStar); | |
1003 | ||
1004 | fPtPion = new TH1F("pionpt","Primary pions candidates pt ",500,0,5); | |
1005 | fPtPion->SetStats(kTRUE); | |
1006 | fPtPion->GetXaxis()->SetTitle("GeV/c"); | |
1007 | fPtPion->GetYaxis()->SetTitle("Entries"); | |
1008 | fOutput->Add(fPtPion); | |
1009 | ||
1010 | fcharmpt = new TH1F("charmpt","MC Charm pt distribution", 10000,0,5000); | |
1011 | fdstarE = new TH1F("dstare", "MC D* energy distribution", 10000,0,5000); | |
1012 | fdstarpt = new TH1F("dstarpt","MC D* pt distribution", 10000,0,5000); | |
1013 | fOutput->Add(fcharmpt); | |
1014 | fOutput->Add(fdstarE); | |
1015 | fOutput->Add(fdstarpt); | |
1016 | ||
1017 | // jet related fistograms | |
1018 | fEjet = new TH1F("ejet", "UA1 algorithm jet energy distribution",1000,0,500); | |
1019 | fPhijet = new TH1F("Phijet","UA1 algorithm jet #phi distribution", 200,-7,7); | |
1020 | fEtaJet = new TH1F("Etajet","UA1 algorithm jet #eta distribution", 200,-7,7); | |
1021 | fOutput->Add(fEjet); | |
1022 | fOutput->Add(fPhijet); | |
1023 | fOutput->Add(fEtaJet); | |
954ac830 | 1024 | |
1025 | fResZ = new TH1F("FragFunc","Fragmentation function ",50,0,1); | |
1026 | fResZBkg = new TH1F("FragFuncBkg","Fragmentation function background",50,0,1); | |
1027 | fOutput->Add(fResZ); | |
1028 | fOutput->Add(fResZBkg); | |
1029 | ||
1030 | fD0ptvsSoftPt = new TH2F("D0piRec","Candidates (background + signal)",100,0,3,100,0,15); | |
1031 | fD0ptvsSoftPtSignal = new TH2F("D0PiSignal","Signal",100,0,3,100,0,15); | |
1032 | fOutput->Add(fD0ptvsSoftPt); | |
1033 | fOutput->Add(fD0ptvsSoftPtSignal); | |
1034 | ||
1035 | return kTRUE; | |
1036 | ||
1037 | } | |
1038 | ||
1039 | //______________________________Phase space cut alternative to PPR _______________________________________ | |
1040 | ||
1041 | Bool_t AliAnalysisTaskSEDStarJets::EvaluateCutOnPiD0pt(AliAODRecoDecayHF2Prong* const vtx, AliAODTrack* const aodTrack) { | |
1042 | ||
1043 | // The soft pion pt and DO pt are strongly correlated. It can be shown that ~ 95% of the signal in constrained | |
1044 | // into a narrow band defined by 10 < D0pt/softPt < 18. This cut can be used with a relaxed CosThetaStar cut | |
1045 | // to reconstruct the D*. | |
1046 | ||
1047 | Double_t softPt = 0; | |
1048 | Double_t d0ptReco = 0; | |
1049 | ||
1050 | softPt = aodTrack->Pt(); | |
1051 | d0ptReco = vtx->Pt(); | |
1052 | fD0ptvsSoftPt->Fill(softPt,d0ptReco); | |
1053 | ||
1054 | if(softPt>0){ | |
1055 | Double_t ratio = d0ptReco/softPt; | |
1056 | if( ratio <=10. || ratio>=18. ) return kFALSE; | |
1057 | } | |
1058 | ||
1059 | fD0ptvsSoftPtSignal->Fill(softPt,d0ptReco); | |
1060 | return kTRUE; | |
1061 | } | |
1062 | ||
1063 | //______________________________ side band background for D*___________________________________ | |
1064 | ||
feaae220 | 1065 | void AliAnalysisTaskSEDStarJets::SideBandBackground(Double_t invM, Double_t invMDStar, Double_t ejet, Double_t dPhi, Int_t nJets){ |
954ac830 | 1066 | |
1067 | // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas | |
1068 | // (expected detector resolution) on the left and right frm the D0 mass. Each band | |
1069 | // has a width of ~5 sigmas. Two band needed for opening angle considerations | |
1070 | ||
1071 | if((invM>=1.763 && invM<=1.811) || (invM>=1.919 && invM<=1.963)){ | |
1072 | ||
feaae220 | 1073 | if(nJets==1)fDiffSideBand->Fill(invMDStar-invM); // M(Kpipi)-M(Kpi) side band background |
954ac830 | 1074 | |
feaae220 | 1075 | if ((invMDStar-invM)<=0.148 && (invMDStar-invM)>=0.142) { |
954ac830 | 1076 | fPhiBkg->Fill(dPhi); |
1077 | ||
1078 | if(dPhi>=-0.5 && dPhi<=0.5){ // evaluate in the near side | |
1079 | ||
1080 | Double_t dStarMomBkg = (fLorentzTrack3 + fLorentzTrack4).P(); | |
1081 | Double_t zFragBkg = (TMath::Cos(dPhi)*dStarMomBkg)/ejet; | |
1082 | fResZBkg->Fill(TMath::Abs(zFragBkg)); | |
1083 | ||
1084 | } | |
1085 | } | |
1086 | } | |
1087 | } |