]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliAnalysisTaskCharmFraction.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
I improved (by a factor 2.5) the speed of the MatchToMC method
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskCharmFraction.cxx
CommitLineData
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
47ClassImp(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//________________________________________________________________________
103AliAnalysisTaskCharmFraction::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//________________________________________________________________________
156void 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//________________________________________________________________________
189void 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//________________________________________________________________________
259void 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//________________________________________________________________________
664void 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}