]>
Commit | Line | Data |
---|---|---|
0df0132a | 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 | // | |
18 | // Class AliAnalysisTaskCharmFraction | |
19 | // AliAnalysisTask for the extraction of the fraction of prompt charm | |
20 | // using the charm hadron impact parameter to the primary vertex | |
21 | // | |
22 | // Author: Andrea Rossi, andrea.rossi@ts.infn.it | |
23 | ///////////////////////////////////////////////////////////// | |
24 | ||
25 | #include <TChain.h> | |
26 | #include <TTree.h> | |
27 | #include <TH1F.h> | |
28 | #include <TH2F.h> | |
29 | #include <TCanvas.h> | |
30 | #include <TDatabasePDG.h> | |
31 | #include <TMath.h> | |
32 | ||
33 | #include "AliAnalysisTask.h" | |
34 | #include "AliAnalysisManager.h" | |
35 | #include "AliAODEvent.h" | |
36 | #include "AliAODInputHandler.h" | |
37 | #include "AliAnalysisVertexingHF.h" | |
38 | #include "AliAODRecoDecayHF2Prong.h" | |
39 | #include "AliAODRecoDecayHF.h" | |
40 | #include "AliAODRecoDecay.h" | |
41 | #include "AliAODTrack.h" | |
42 | #include "AliAODVertex.h" | |
43 | #include "AliAODMCParticle.h" | |
44 | #include "AliAODMCHeader.h" | |
45 | #include "AliAnalysisTaskCharmFraction.h" | |
46 | ||
47 | ClassImp(AliAnalysisTaskCharmFraction) | |
48 | ||
49 | //________________________________________________________________________ | |
50 | AliAnalysisTaskCharmFraction::AliAnalysisTaskCharmFraction(const char *name) | |
51 | : AliAnalysisTask(name, ""), fAOD(0), fVHF(0), | |
52 | fArrayD0toKpi(0), | |
53 | fArrayMC(0), | |
54 | fAODmcHeader(0), | |
55 | fhMass(0), | |
56 | fhMassTrue(0), | |
57 | fhCPtaVSd0d0(0), | |
58 | fhd0xd0(0), | |
59 | fhCPta(0), | |
60 | fhSecVtxZ(0), | |
61 | fhSecVtxX(0), | |
62 | fhSecVtxY(0), | |
63 | fhSecVtxXY(0), | |
64 | fhSecVtxPhi(0), | |
65 | fhd0D0(0), | |
66 | fhd0D0pt(0x0), | |
67 | fhd0D0VtxTrue(0), | |
68 | fhd0D0VtxTruept(0x0), | |
69 | fhMCd0D0(0), | |
70 | fhMCd0D0pt(0x0), | |
71 | fnbins(), | |
72 | fD0usecuts(kTRUE), | |
73 | fcheckMC(kFALSE), | |
74 | fcheckMCD0(kFALSE), | |
75 | fcheckMC2prongs(kFALSE), | |
76 | fcheckMCprompt(kFALSE), | |
77 | fcheckMCfromB(kFALSE), | |
78 | fSkipD0star(kFALSE), | |
79 | fStudyd0fromBTrue(kTRUE), | |
80 | fStudyPureBackground(kFALSE), | |
81 | fSideBands(0) | |
82 | { | |
83 | // Constructor | |
84 | //Side bands selection (see cxx): 999 -> no inv mass selection at all, | |
85 | // >0 -> both D0 and D0bar hypothesis are required to be kFALSE, | |
86 | // <0 -> only 1 inv mass hypothesis is required to be KFALSE | |
87 | ||
88 | SetNPtBins(); | |
89 | // Define input and output slots here | |
90 | // Input slot #0 works with a TChain | |
91 | DefineInput(0, TChain::Class()); | |
92 | // DefineInput(1, TChain::Class()); | |
93 | // Output slot #0 writes into a TH1 container | |
94 | DefineOutput(0, TH2F::Class());// hCPtaVsd0d0 | |
95 | DefineOutput(1, TH2F::Class());// hVtxXY | |
96 | for(Int_t j=2;j<3*fnbins+11;j++){ // hCPta, hd0xd0, hVtZ, hVtX hVtY hVtPhi + nptbin hd0D0 +hd0D0 integreated +nptbin hd0VtxTrueD0 +hd0VtxTrueD0 integrated +nptbin hMCd0D0 +hMCd0D0 integrated | |
97 | DefineOutput(j, TH1F::Class()); | |
98 | } | |
99 | // DefineOutput(fnbins+8, TH1F::Class());// hd0D0 pt integrated | |
100 | } | |
101 | ||
102 | //________________________________________________________________________ | |
103 | AliAnalysisTaskCharmFraction::AliAnalysisTaskCharmFraction(const char *name,Int_t nptbins) | |
104 | : AliAnalysisTask(name, ""), fAOD(0), fVHF(0), | |
105 | fArrayD0toKpi(0), | |
106 | fArrayMC(0), | |
107 | fAODmcHeader(0), | |
108 | fhMass(0), | |
109 | fhMassTrue(0), | |
110 | fhCPtaVSd0d0(0), | |
111 | fhd0xd0(0), | |
112 | fhCPta(0), | |
113 | fhSecVtxZ(0), | |
114 | fhSecVtxX(0), | |
115 | fhSecVtxY(0), | |
116 | fhSecVtxXY(0), | |
117 | fhSecVtxPhi(0), | |
118 | fhd0D0(0), | |
119 | fhd0D0pt(0x0), | |
120 | fhd0D0VtxTrue(0), | |
121 | fhd0D0VtxTruept(0x0), | |
122 | fhMCd0D0(0), | |
123 | fhMCd0D0pt(0x0), | |
124 | fnbins(), | |
125 | fD0usecuts(kTRUE), | |
126 | fcheckMC(kFALSE), | |
127 | fcheckMCD0(kFALSE), | |
128 | fcheckMC2prongs(kFALSE), | |
129 | fcheckMCprompt(kFALSE), | |
130 | fcheckMCfromB(kFALSE), | |
131 | fSkipD0star(kFALSE), | |
132 | fStudyd0fromBTrue(kTRUE), | |
133 | fStudyPureBackground(kFALSE), | |
134 | fSideBands(0) | |
135 | { | |
136 | // Constructor | |
137 | //Side bands selection (see cxx): 999 -> no inv mass selection at all, | |
138 | // >0 -> both D0 and D0bar hypothesis are required to be kFALSE, | |
139 | // <0 -> only 1 inv mass hypothesis is required to be KFALSE | |
140 | SetNPtBins(nptbins); | |
141 | // Define input and output slots here | |
142 | // Input slot #0 works with a TChain | |
143 | DefineInput(0, TChain::Class()); | |
144 | // DefineInput(1, TChain::Class()); | |
145 | // Output slot #0 writes into a TH1 container | |
146 | ||
147 | DefineOutput(0, TH2F::Class());// hCPtaVsd0d0 | |
148 | DefineOutput(1, TH2F::Class());// hVtxXY | |
149 | for(Int_t j=2;j<3*fnbins+13;j++){ // hCPta, hd0xd0, hVtZ, hVtX hVtY hVtPhi + nptbin hd0D0 +hd0D0 integreated +nptbin hd0VtxTrueD0 +hd0VtxTrueD0 integrated +nptbin hMCd0D0 +hMCd0D0 integrated | |
150 | DefineOutput(j, TH1F::Class()); | |
151 | } | |
152 | // DefineOutput(fnbins+8, TH1F::Class());// hd0D0 pt integrated | |
153 | } | |
154 | ||
155 | //________________________________________________________________________ | |
156 | void AliAnalysisTaskCharmFraction::ConnectInputData(Option_t *) | |
157 | { | |
158 | // Connect ESD or AOD here | |
159 | // Called once | |
160 | ||
161 | TTree* tree = dynamic_cast<TTree*> (GetInputData(0)); | |
162 | // TTree* treeFriend = dynamic_cast<TTree*> (GetInputData(1)); | |
163 | if (!tree) { | |
164 | Printf("ERROR: Could not read chain from input slot 0"); | |
165 | } else { | |
166 | // tree->AddFriend(treeFriend); | |
167 | ||
168 | // Disable all branches and enable only the needed ones | |
169 | // The next two lines are different when data produced as AliESDEvent is read | |
170 | ||
171 | AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); | |
172 | fAOD=new AliAODEvent(); | |
173 | fAOD->ReadFromTree(tree); | |
174 | ||
175 | fArrayD0toKpi = (TClonesArray*)fAOD->GetList()->FindObject("D0toKpi"); | |
176 | fArrayMC = (TClonesArray*)fAOD->GetList()->FindObject("mcparticles"); | |
177 | ||
178 | fAODmcHeader=(AliAODMCHeader*)fAOD->GetList()->FindObject("mcHeader"); | |
179 | ||
180 | if (!aodH) { | |
181 | Printf("ERROR: Could not get AODInputHandler"); | |
182 | } else | |
183 | fAOD = aodH->GetEvent(); | |
184 | ||
185 | } | |
186 | } | |
187 | ||
188 | //________________________________________________________________________ | |
189 | void AliAnalysisTaskCharmFraction::CreateOutputObjects() | |
190 | { | |
191 | // Create histograms | |
192 | // Called once | |
193 | fhd0D0 = new TH1F("hd0D0","D^{0} impact par. plot (All momenta)",1000,-1000.,1000.); | |
194 | fhd0D0->SetXTitle("Impact parameter [#mum]"); | |
195 | fhd0D0->SetYTitle("Entries"); | |
196 | ||
197 | fhd0D0VtxTrue = new TH1F("hd0D0VtxTrue","D^{0} impact par. w.r.t. True Vtx (All momenta)",1000,-1000.,1000.); | |
198 | fhd0D0VtxTrue->SetXTitle("Impact parameter [#mum]"); | |
199 | fhd0D0VtxTrue->SetYTitle("Entries"); | |
200 | ||
201 | fhMCd0D0 = new TH1F("hMCd0D0","D^{0} impact par. plot (All momenta)",1000,-1000.,1000.); | |
202 | fhMCd0D0->SetXTitle("MC Impact parameter [#mum]"); | |
203 | fhMCd0D0->SetYTitle("Entries"); | |
204 | ||
205 | fhCPtaVSd0d0=new TH2F("hCPtaVSd0d0","hCPtaVSd0d0",1000,-100000.,100000.,100,0.,1.); | |
206 | fhSecVtxZ=new TH1F("hSecVtxZ","hSecVtxZ",1000,-8.,8.); | |
207 | fhSecVtxX=new TH1F("hSecVtxX","hSecVtxX",1000,-3000.,3000.); | |
208 | fhSecVtxY=new TH1F("hSecVtxY","hSecVtxY",1000,-3000.,3000.); | |
209 | fhSecVtxXY=new TH2F("hSecVtxXY","hSecVtxXY",1000,-3000.,3000.,1000,-3000.,3000.); | |
210 | fhSecVtxPhi=new TH1F("hSecVtxPhi","hSecVtxPhi",180,-180.1,180.1); | |
211 | fhCPta=new TH1F("hCPta","hCPta",100,0.,1.); | |
212 | fhd0xd0=new TH1F("hd0xd0","hd0xd0",1000,-100000.,100000.); | |
213 | fhMassTrue=new TH1F("hMassTrue","D^{0} MC inv. Mass (All momenta)",600,1.600,2.200); | |
214 | fhMass=new TH1F("hMass","D^{0} inv. Mass (All momenta)",600,1.600,2.200); | |
215 | Double_t standardptbins[14]={0.,0.5,1.,2.,3.,5.,7.,10.,15.,20.,30.,40.,50.,100.}; | |
216 | ||
217 | fhd0D0pt=new TH1F*[fnbins]; | |
218 | fhMCd0D0pt=new TH1F*[fnbins]; | |
219 | fhd0D0VtxTruept=new TH1F*[fnbins]; | |
220 | TString namehist="hd0D0pt"; | |
221 | TString titlehist="D^{0} impact par. plot "; | |
222 | TString strnamept,strtitlept; | |
223 | for(Int_t i=0;i<fnbins;i++){ | |
224 | strnamept=namehist; | |
225 | strnamept+=i;// strnamept+=standardptbins[i+1]; | |
226 | // strnamept.Append("_"); | |
227 | //strnamept+=standardptbins[i+1]; | |
228 | ||
229 | strtitlept=namehist; | |
230 | strtitlept+=standardptbins[i]; | |
231 | strtitlept.Append("<= pt <"); | |
232 | strtitlept+=standardptbins[i+1]; | |
233 | strtitlept.Append(" [GeV/c]"); | |
234 | ||
235 | fhd0D0pt[i] = new TH1F(strnamept.Data(),strtitlept.Data(),1000,-1000.,1000.); | |
236 | fhd0D0pt[i]->SetXTitle("Impact parameter [#mum] "); | |
237 | fhd0D0pt[i]->SetYTitle("Entries"); | |
238 | ||
239 | strnamept.ReplaceAll("hd0D0","hMCd0D0"); | |
240 | fhMCd0D0pt[i] = new TH1F(strnamept.Data(),strtitlept.Data(),1000,-1000.,1000.); | |
241 | fhMCd0D0pt[i]->SetXTitle("MC Impact parameter [#mum] "); | |
242 | fhMCd0D0pt[i]->SetYTitle("Entries"); | |
243 | ||
244 | strnamept.ReplaceAll("hMCd0D0","hd0D0VtxTrue"); | |
245 | fhd0D0VtxTruept[i] = new TH1F(strnamept.Data(),strtitlept.Data(),1000,-1000.,1000.); | |
246 | fhd0D0VtxTruept[i]->SetXTitle("Impact parameter w.r.t. True Vtx [#mum] "); | |
247 | fhd0D0VtxTruept[i]->SetYTitle("Entries"); | |
248 | ||
249 | } | |
250 | ||
251 | ||
252 | /* fHistPt = new TH1F("fHistPt", "P_{T} distribution", 15, 0.1, 3.1); | |
253 | fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)"); | |
254 | fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); | |
255 | fHistPt->SetMarkerStyle(kFullCircle);*/ | |
256 | } | |
257 | ||
258 | //________________________________________________________________________ | |
259 | void AliAnalysisTaskCharmFraction::Exec(Option_t *) | |
260 | { | |
261 | // Main loop | |
262 | // Called for each event | |
263 | ||
264 | if (!fAOD) { | |
265 | Printf("ERROR: fAOD not available"); | |
266 | return; | |
267 | } | |
268 | ||
269 | //Printf("There are %d tracks in this event", fAOD->GetNumberOfTracks()); | |
270 | // Int_t nTotHF=0,nTotDstar=0,nTot3Prong=0; | |
271 | Int_t nTotD0toKpi=0; | |
272 | // const Int_t nptbins=10; | |
273 | Double_t standardptbins[14]={0.,0.5,1.,2.,3.,5.,7.,10.,15.,20.,30.,40.,50.,100.}; | |
274 | Double_t pD[3],xD[3],pXtrTrue[2],pYtrTrue[2],pZtrTrue[2],xtr1[3],xtr2[3]; | |
275 | Double_t cutsD0[9]= | |
276 | // cutsD0[0] = inv. mass half width [GeV] | |
277 | // cutsD0[1] = dca [cm] | |
278 | // cutsD0[2] = cosThetaStar | |
279 | // cutsD0[3] = pTK [GeV/c] | |
280 | // cutsD0[4] = pTPi [GeV/c] | |
281 | // cutsD0[5] = d0K [cm] upper limit! | |
282 | // cutsD0[6] = d0Pi [cm] upper limit! | |
283 | // cutsD0[7] = d0d0 [cm^2] | |
284 | // cutsD0[8] = cosThetaPoint | |
285 | {1000., | |
286 | 100000., | |
287 | 1.1, | |
288 | 0., | |
289 | 0., | |
290 | 100000., | |
291 | 100000., | |
292 | 100000000., | |
293 | -1.1}; | |
294 | if(fD0usecuts){ | |
295 | cutsD0[0]=0.033; | |
296 | cutsD0[1]=0.03; | |
297 | cutsD0[2]=0.8; | |
298 | cutsD0[3]=0.; | |
299 | cutsD0[4]=0.; | |
300 | cutsD0[5]=100000.; | |
301 | cutsD0[6]=100000.; | |
302 | cutsD0[7]=-0.00005; | |
303 | cutsD0[8]=0.8; | |
304 | } | |
305 | const Double_t fixedInvMassCut=cutsD0[0]; | |
306 | AliAODVertex *vtx1=0; | |
307 | // primary vertex | |
308 | vtx1 = (AliAODVertex*)fAOD->GetPrimaryVertex(); | |
309 | Double_t vtxTrue[3]; | |
310 | fAODmcHeader->GetVertex(vtxTrue); | |
311 | // vtx1->Print(); | |
312 | AliAODRecoDecayHF *aodDMC=0x0; | |
313 | Bool_t nomum=kFALSE,background=kFALSE; | |
314 | AliAODMCParticle *mum1=0x0,*mum2=0x0; | |
315 | AliAODMCParticle *b1=0x0,*b2=0x0; | |
316 | AliAODMCParticle *sonpart=0x0,*grandmoth1=0x0; | |
317 | ||
318 | // make trkIDtoEntry register (temporary) | |
319 | Int_t trkIDtoEntry[100000]; | |
320 | for(Int_t it=0;it<fAOD->GetNumberOfTracks();it++) { | |
321 | AliAODTrack *track = fAOD->GetTrack(it); | |
322 | if(track->GetID()<0) { | |
323 | //printf("Track ID <0, id= %d\n",track->GetID()); | |
324 | return; | |
325 | } | |
326 | trkIDtoEntry[track->GetID()]=it; | |
327 | } | |
328 | ||
329 | ||
330 | // loop over D0->Kpi candidates | |
331 | Int_t nD0toKpi = fArrayD0toKpi->GetEntriesFast(); | |
332 | nTotD0toKpi += nD0toKpi; | |
333 | // cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl; | |
334 | ||
335 | for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) { | |
336 | if(aodDMC!=0x0)delete aodDMC; | |
337 | mum1=0x0; | |
338 | mum2=0x0; | |
339 | b1=0x0; | |
340 | b2=0x0; | |
341 | sonpart=0x0; | |
342 | grandmoth1=0x0; | |
343 | ||
344 | AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)fArrayD0toKpi->UncheckedAt(iD0toKpi); | |
345 | Bool_t unsetvtx=kFALSE; | |
346 | if(!d->GetOwnPrimaryVtx()) { | |
347 | d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables | |
348 | unsetvtx=kTRUE; | |
349 | } | |
350 | ||
351 | // d->SetOwnPrimaryVtx(vtx1); | |
352 | //unsetvtx=kTRUE; | |
353 | ||
354 | Int_t okD0=0,okD0bar=0; | |
438f219a | 355 | Double_t invMassD0,invMassD0bar,cutmassvalue=-1.; |
0df0132a | 356 | Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass(); |
357 | nomum=kFALSE; | |
358 | background=kFALSE; | |
359 | ||
360 | cutsD0[0]=fixedInvMassCut; | |
361 | if(fSideBands){ | |
362 | cutmassvalue=cutsD0[0]; | |
363 | if(fSideBands==999.)cutsD0[0]=1000.; | |
364 | else cutsD0[0]=2.*TMath::Abs(fSideBands)*cutsD0[0]; | |
365 | } | |
366 | if(d->SelectD0(cutsD0,okD0,okD0bar)) { | |
367 | if(fSideBands){ | |
368 | if(fSideBands==999.)cutmassvalue=-1.; | |
369 | d->InvMassD0(invMassD0,invMassD0bar); | |
370 | ||
371 | if(TMath::Abs(invMassD0-mD0PDG) > TMath::Abs(fSideBands)*cutmassvalue) okD0 = okD0 && kTRUE; | |
372 | else okD0=kFALSE; | |
373 | if(TMath::Abs(invMassD0bar-mD0PDG) > TMath::Abs(fSideBands)*cutmassvalue) okD0bar = okD0bar && kTRUE; // SIDE BANDS SELECTION | |
374 | else okD0bar=kFALSE; | |
375 | ||
376 | if(fSideBands<0.){ | |
377 | if(!okD0 && !okD0bar) continue; | |
378 | } | |
379 | else if (!okD0 || !okD0bar)continue; | |
380 | } | |
381 | ||
382 | ||
383 | // get daughter AOD tracks | |
384 | AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0); | |
385 | AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1); | |
386 | ||
387 | //Check if it's a True D0 | |
388 | if(fcheckMC){ | |
389 | while(!background){ | |
390 | ||
391 | if(trk0->GetLabel()==-1||trk1->GetLabel()==-1){ | |
392 | //fake tracks!! | |
393 | background=kTRUE; //FAKE TRACK | |
394 | break; | |
395 | } | |
396 | // printf("Before entering the MC checks \n"); | |
397 | ||
398 | b1=(AliAODMCParticle*)fArrayMC->At(trk0->GetLabel()); | |
399 | b2=(AliAODMCParticle*)fArrayMC->At(trk1->GetLabel()); | |
400 | ||
401 | ||
402 | if(b1->GetMother()==-1||b2->GetMother()==-1){ | |
403 | nomum=kTRUE; //Tracks with no mother | |
404 | //printf("Inside problem mothering \n");// FAKE DECAY VERTEX | |
405 | background=kTRUE; | |
406 | break; | |
407 | // if(!fStudyPureBackground)continue; | |
408 | } | |
409 | ||
410 | if(b1->GetMother()!=b2->GetMother()){ | |
411 | //Check the label of the mother is the same | |
412 | background=kTRUE; // NOT SAME MOTHER | |
413 | break; | |
414 | } | |
415 | // printf("After first mothering \n"); | |
416 | mum1=(AliAODMCParticle*)fArrayMC->At(b1->GetMother()); | |
417 | mum2=(AliAODMCParticle*)fArrayMC->At(b2->GetMother()); | |
418 | // if((mum1->GetPdgCode()!=mum2->GetPdgCode()))continue; //Check the mother is the same particle | |
419 | // printf("After second mothering \n"); | |
420 | // printf("Particle codes: tr1: %d, tr2: %d, mum1: %d, mum 2: %d \n",b1->GetPdgCode(),b2->GetPdgCode(),mum1->GetPdgCode(),mum2->GetPdgCode()); | |
421 | if(fcheckMCD0)if(TMath::Abs(mum1->GetPdgCode())!=421){ | |
422 | //^fStudyPureBackground)continue;//NPRONG | |
423 | background=kTRUE;// NOT A D0 | |
424 | break; | |
425 | } | |
426 | if(fcheckMC2prongs)if(mum1->GetDaughter(1)-mum1->GetDaughter(0)+1!=2){ | |
427 | //^fStudyPureBackground)continue;// NPRONG | |
428 | background=kTRUE; // NOT A 2 PRONG DECAY | |
429 | break; | |
430 | } | |
431 | ||
432 | if(mum1->GetMother()==-1){ | |
433 | background=kTRUE; // A particle coming from nothing | |
434 | break; | |
435 | } | |
436 | ||
437 | // ^fStudyPureBackground)continue; | |
438 | // printf("After second mothering \n"); | |
439 | grandmoth1=(AliAODMCParticle*)fArrayMC->At(mum1->GetMother()); | |
440 | if(fSkipD0star)if(TMath::Abs(grandmoth1->GetPdgCode())==413||TMath::Abs(grandmoth1->GetPdgCode())==423){ | |
441 | //^fStudyPureBackground)continue; | |
442 | background=kTRUE;// D0 COMING FROM A D0* | |
443 | break; | |
444 | } | |
445 | if(fcheckMCD0){//THIS CHECK IS NEEDE TO AVOID POSSIBLE FAILURE IN THE SECOND WHILE | |
446 | while(TMath::Abs(grandmoth1->GetPdgCode())!=4&&TMath::Abs(grandmoth1->GetPdgCode())!=5){ | |
447 | if(grandmoth1->GetMother()==-1){ | |
448 | // printf("Inside grandmothering \n"); | |
449 | //printf("mother=-1, pdgcode: %d \n",grandmoth1->GetPdgCode()); | |
450 | Int_t son=grandmoth1->GetDaughter(0); | |
451 | sonpart=(AliAODMCParticle*)fArrayMC->At(son); | |
452 | while(TMath::Abs(sonpart->GetPdgCode())!=421){ | |
453 | //printf("mother=-1, pdgcode: %d \n",sonpart->GetPdgCode()); | |
454 | son++; | |
455 | sonpart=(AliAODMCParticle*)fArrayMC->At(son); | |
456 | } | |
457 | break; | |
458 | } | |
459 | grandmoth1=(AliAODMCParticle*)fArrayMC->At(grandmoth1->GetMother()); | |
460 | } | |
461 | } | |
462 | if(fcheckMCprompt)if(TMath::Abs(grandmoth1->GetPdgCode())!=4){ | |
463 | //^fStudyPureBackground)continue; | |
464 | background=kTRUE; | |
465 | break; | |
466 | } | |
467 | if(fcheckMCfromB){ | |
468 | if(TMath::Abs(grandmoth1->GetPdgCode())!=5){ | |
469 | //^fStudyPureBackground)continue; | |
470 | background=kTRUE; | |
471 | break; | |
472 | } | |
473 | } | |
474 | //else if(!fStudyPureBackground)continue; | |
475 | break; | |
476 | } | |
477 | if(background^fStudyPureBackground)continue; | |
478 | if(fStudyd0fromBTrue){ | |
479 | if(b1!=0x0&&b2!=0x0){ | |
480 | Int_t charge[2]={0,0}; | |
481 | if(b1->Charge()==-1)charge[0]=1; | |
482 | else { | |
483 | if(b2->Charge()==-1){ | |
484 | //printf("Same charges for prongs \n"); | |
485 | continue; | |
486 | } | |
487 | charge[1]=1; | |
488 | } | |
489 | /* if(!pXtrTrue[charge[0]]=b1->Px()){printf("!first MCtrack:Get momentum X failed \n");continue;} | |
490 | if(!b1->Py(pYtrTrue[charge[0]])){printf("!first MCtrack:Get momentum Y failed \n");continue;} | |
491 | if(!b1->Pz(pZtrTrue[charge[0]])){printf("!first MCtrack:Get momentum Z failed \n");continue;} | |
492 | */ | |
493 | ||
494 | pXtrTrue[charge[0]]=b1->Px(); | |
495 | pYtrTrue[charge[0]]=b1->Py(); | |
496 | pZtrTrue[charge[0]]=b1->Pz(); | |
497 | if(!b1->XvYvZv(xtr1)){ | |
498 | //printf("!first MCtrack:Get position failed \n"); | |
499 | continue; | |
500 | } | |
501 | /*if(!b2->Px(pXtrTrue[charge[1]])){printf("!second MCtrack:Get momentum X failed \n");continue;} | |
502 | if(!b2->Py(pYtrTrue[charge[1]])){printf("!second MCtrack:Get momentum Y failed \n");continue;} | |
503 | if(!b2->Pz(pZtrTrue[charge[1]])){printf("!second MCtrack:Get momentum Z failed \n");continue;}*/ | |
504 | pXtrTrue[charge[1]]=b2->Px(); | |
505 | pYtrTrue[charge[1]]=b2->Py(); | |
506 | pZtrTrue[charge[1]]=b2->Pz(); | |
507 | if(!nomum){ | |
508 | if(!(mum1==0x0||mum2==0x0)){ | |
509 | if(!mum1->PxPyPz(pD)){ | |
510 | //printf("!D from B:Get momentum failed \n"); | |
511 | continue; | |
512 | } | |
513 | if(!mum1->XvYvZv(xD)){ | |
514 | //printf("!D from B:Get position failed \n"); | |
515 | continue; | |
516 | } | |
517 | if(pXtrTrue[0]+pXtrTrue[1]!=pD[0]){ | |
518 | //printf("Warning! MCinfo: Pt of the mother doesn't match the sum of the daughters pt: X component %f + %f = %f Vs %f \n",pXtrTrue[0],pXtrTrue[1],pXtrTrue[0]+pXtrTrue[1],pD[0]); | |
519 | //printf("Warning! MCinfo: Pt of the mother doesn't match the sum of the daughters pt: Y componente %f + %f = %f Vs %f \n",pYtrTrue[0],pYtrTrue[1],pYtrTrue[0]+pYtrTrue[1],pD[1]); | |
520 | //printf("Warning!: Pt comparison from the sum of the daughters:%f Vs mother pT: %f \n",TMath::Sqrt((pXtrTrue[0]+pXtrTrue[1])*(pXtrTrue[0]+pXtrTrue[1])+(pYtrTrue[0]+pYtrTrue[1])*(pYtrTrue[0]+pYtrTrue[1])),TMath::Sqrt(pD[0]*pD[0]+pD[1]*pD[1])); | |
521 | // return; | |
522 | } | |
523 | ||
524 | // if(!b2->PxPyPz(ptr2True)){printf("!second MCtrack:Get momentum failed \n");continue;} | |
525 | if(!b2->XvYvZv(xtr2)){ | |
526 | //printf("!second MCtrack:Get position failed \n"); | |
527 | continue; | |
528 | } | |
529 | Double_t d0dummy[2]={0.,0.};//TEMPORARY : dummy d0 for AliAODRecoDeay constructor | |
530 | aodDMC=new AliAODRecoDecayHF(vtxTrue,xD,2,0,pXtrTrue,pYtrTrue,pZtrTrue,d0dummy); | |
531 | fhMassTrue->Fill(mum1->GetCalcMass()); | |
532 | } | |
533 | //else printf("Warning! The candidate comes from two track with a different mother, \n -> the true mother d0 doesn't exist ! -> There will be different entries between observed and MC d0distributions \n"); | |
534 | } | |
535 | } | |
536 | } | |
537 | //printf("End of MC checks \n"); | |
538 | } | |
539 | ||
540 | //printf("A real D0 was found!! \n"); | |
541 | ||
542 | ||
543 | //cout<<1e8*d->Prodd0d0()<<endl; | |
544 | /* if(fSideBands){ // The integral of the Hists will in general differ | |
545 | // from the number of candidates which passed the selection | |
546 | if(okD0)fhMass->Fill(d->InvMassD0(),1.); | |
547 | if(okD0bar)fhMass->Fill(d->InvMassD0bar(),1.); | |
548 | } | |
549 | else { // The integral of the Hists will in general differ | |
550 | // from the number of candidates which passed the selection | |
551 | if(okD0)fhMass->Fill(d->InvMassD0(),1.); | |
552 | if(okD0bar)fhMass->Fill(d->InvMassD0bar(),1.); | |
553 | } | |
554 | */ | |
555 | ||
556 | if(okD0)fhMass->Fill(d->InvMassD0(),1.); | |
557 | if(okD0bar)fhMass->Fill(d->InvMassD0bar(),1.); | |
558 | ||
559 | fhCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle()); | |
560 | fhSecVtxZ->Fill(d->GetSecVtxZ()); | |
561 | fhSecVtxX->Fill(d->GetSecVtxX()*10000.); | |
562 | fhSecVtxY->Fill(d->GetSecVtxY()*10000.); | |
563 | fhSecVtxXY->Fill(d->GetSecVtxX()*10000.,d->GetSecVtxY()*10000.); | |
564 | fhSecVtxPhi->Fill(TMath::ATan2(d->GetSecVtxY()*10000.,d->GetSecVtxX()*10000.)*TMath::RadToDeg()); | |
565 | fhd0xd0->Fill(1e8*d->Prodd0d0()); | |
566 | fhCPta->Fill(d->CosPointingAngle()); | |
567 | ||
568 | ||
569 | //cout<<d->GetSecVtxX() <<endl; | |
570 | ||
571 | if(!trk0 || !trk1) { | |
572 | /* if(d->GetProngID(0)<0||d->GetProngID(1)<0){ | |
573 | printf("Wrong ProngID,%d or %d\n",d->GetProngID(0),d->GetProngID(1));continue; | |
574 | } | |
575 | */ | |
576 | trk0=fAOD->GetTrack(trkIDtoEntry[d->GetProngID(0)]); | |
577 | trk1=fAOD->GetTrack(trkIDtoEntry[d->GetProngID(1)]); | |
578 | // cout<<"references to standard AOD not available"<<endl; | |
579 | } | |
580 | // cout<<"pt of positive track: "<<trk0->Pt()<<endl; | |
581 | ||
582 | // make a AliExternalTrackParam from the D0 | |
583 | // and calculate impact parameters to primary vertex | |
584 | ||
585 | ||
586 | // AliExternalTrackParam trackD0(d); | |
587 | //cout<<"pt of D0: "<<d->Pt()<<" (AliAODRecoDecay); "<<trackD0.Pt()<<" (track)"<<endl; | |
588 | //trackD0.Print(); | |
589 | Double_t ptD0=d->Pt(); | |
590 | for(Int_t j=0;j<fnbins;j++){ | |
591 | if(ptD0<standardptbins[j+1]){ | |
592 | fhd0D0VtxTruept[j]->Fill(d->AliAODRecoDecay::ImpParXY(vtxTrue)*10000.); | |
593 | fhd0D0pt[j]->Fill(d->ImpParXY()*10000.); | |
594 | // printf("HERE in tha class \n"); | |
595 | if(aodDMC!=0x0)fhMCd0D0pt[j]->Fill(aodDMC->ImpParXY()*10000.); | |
596 | // printf("HERE2 in tha class \n"); | |
597 | break; | |
598 | } | |
599 | } | |
600 | fhd0D0->Fill(d->ImpParXY()*10000.); | |
601 | fhd0D0VtxTrue->Fill(d->AliAODRecoDecay::ImpParXY(vtxTrue)*10000.); | |
602 | if(aodDMC!=0x0)fhMCd0D0->Fill(aodDMC->ImpParXY()*10000.); | |
603 | // if((1.864-3.*0.011.<d->InvMassD0()&&d->InvMassD0()<1.864+3.*0.011)||(1.864-3.*0.011<d->InvMassD0bar()&&d->InvMassD0bar()<1.864+3.*0.011)){ | |
604 | if(aodDMC!=0x0){ | |
605 | delete aodDMC; | |
606 | aodDMC=0x0; | |
607 | } | |
608 | mum1=0x0; | |
609 | mum2=0x0; | |
610 | b1=0x0; | |
611 | b2=0x0; | |
612 | sonpart=0x0; | |
613 | grandmoth1=0x0; | |
614 | ||
615 | // } | |
616 | // Double_t dz[2],covdz[3]; | |
617 | //trackD0.PropagateToDCA(vtx1,fAOD->GetMagneticField(),1000.,dz,covdz); | |
618 | // cout<<"D0 impact parameter rphi: "<<dz[0]<<" +- "<<TMath::Sqrt(covdz[0])<<endl; | |
619 | } | |
620 | // printf("Before end of the event \n"); | |
621 | if(unsetvtx) d->UnsetOwnPrimaryVtx(); | |
622 | ||
623 | } | |
624 | ||
625 | ||
626 | ||
627 | // Post output data. | |
628 | ||
629 | PostData(0,fhCPtaVSd0d0); | |
630 | PostData(1,fhSecVtxXY); | |
631 | PostData(2,fhd0xd0); | |
632 | PostData(3,fhCPta); | |
633 | PostData(4,fhSecVtxZ); | |
634 | PostData(5,fhSecVtxX); | |
635 | PostData(6,fhSecVtxY); | |
636 | PostData(7,fhSecVtxPhi); | |
637 | ||
638 | ||
639 | for(Int_t j=0;j<fnbins;j++){ | |
640 | ||
641 | PostData(j+8, fhd0D0pt[j]); | |
642 | } | |
643 | ||
644 | PostData(fnbins+8, fhd0D0); | |
645 | ||
646 | for(Int_t j=0;j<fnbins;j++){ | |
647 | ||
648 | PostData(j+fnbins+9,fhMCd0D0pt[j]); | |
649 | } | |
650 | ||
651 | PostData(2*fnbins+9, fhMCd0D0); | |
652 | ||
653 | for(Int_t j=0;j<fnbins;j++){ | |
654 | ||
655 | PostData(j+2*fnbins+10,fhd0D0VtxTruept[j]); | |
656 | } | |
657 | ||
658 | PostData(3*fnbins+10,fhd0D0VtxTrue); | |
659 | PostData(3*fnbins+11,fhMassTrue); | |
660 | PostData(3*fnbins+12,fhMass); | |
661 | } | |
662 | ||
663 | //________________________________________________________________________ | |
664 | void AliAnalysisTaskCharmFraction::Terminate(Option_t *) | |
665 | { | |
666 | // Draw result to the screen | |
667 | // Called once at the end of the query | |
668 | ||
669 | fhd0D0 = dynamic_cast<TH1F*> (GetOutputData(fnbins+8)); | |
670 | if (!fhd0D0) { | |
671 | Printf("ERROR: fhd0D0 not available"); | |
672 | return; | |
673 | } | |
674 | ||
675 | TCanvas *c1 = new TCanvas("AliAnalysisTaskCharmFraction","D0 impact par",10,10,510,510); | |
676 | c1->cd(1)->SetLogy(); | |
677 | fhd0D0->DrawCopy("E"); | |
678 | ||
679 | TCanvas *cd0D0pt=new TCanvas("cd0D0pt","cd0D0pt",0,0,800,700); | |
680 | cd0D0pt->Divide(fnbins/2+fnbins%2,2); | |
681 | for(Int_t j=0;j<fnbins;j++){ | |
682 | fhd0D0pt[j] = dynamic_cast<TH1F*> (GetOutputData(8+j)); | |
683 | if (!fhd0D0pt[j]) { | |
684 | Printf("ERROR: fhd0D0pt not available"); | |
685 | return; | |
686 | } | |
687 | cd0D0pt->cd(j+1); | |
688 | fhd0D0pt[j]->Draw(); | |
689 | } | |
690 | ||
691 | } |