]>
Commit | Line | Data |
---|---|---|
3a219f60 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-2008, 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 | // AliAnalysisTaskSE for the comparison of heavy flavor | |
19 | // decay candidates with the MC truth. | |
20 | // | |
21 | // Author: A.Dainese, andrea.dainese@lnl.infn.it | |
22 | ///////////////////////////////////////////////////////////// | |
23 | ||
24 | #include <TClonesArray.h> | |
25 | #include <TNtuple.h> | |
26 | #include <TList.h> | |
27 | #include <TH1F.h> | |
28 | ||
29 | #include "AliAODEvent.h" | |
30 | #include "AliAODVertex.h" | |
607719af | 31 | #include "AliESDtrack.h" |
3a219f60 | 32 | #include "AliAODTrack.h" |
33 | #include "AliAODMCHeader.h" | |
34 | #include "AliAODMCParticle.h" | |
35 | #include "AliAODRecoDecayHF2Prong.h" | |
b3999999 | 36 | #include "AliAODRecoCascadeHF.h" |
2ae9c7cb | 37 | #include "AliAnalysisVertexingHF.h" |
3a219f60 | 38 | #include "AliAnalysisTaskSE.h" |
39 | #include "AliAnalysisTaskSECompareHF.h" | |
40 | ||
41 | ClassImp(AliAnalysisTaskSECompareHF) | |
42 | ||
43 | ||
44 | //________________________________________________________________________ | |
45 | AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(): | |
46 | AliAnalysisTaskSE(), | |
47 | fOutput(0), | |
2a041947 | 48 | fNtupleCmp(0), |
12bfc069 | 49 | fHistMass(0), |
a273622f | 50 | fHistNEvents(0), |
12bfc069 | 51 | fVHF(0) |
3a219f60 | 52 | { |
53 | // Default constructor | |
3a219f60 | 54 | |
9d6d35b0 | 55 | // NO DefineOutput() HERE (ONLY IN STANDARD CONSTRUCTOR) |
3a219f60 | 56 | } |
57 | ||
58 | //________________________________________________________________________ | |
59 | AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name): | |
60 | AliAnalysisTaskSE(name), | |
61 | fOutput(0), | |
2a041947 | 62 | fNtupleCmp(0), |
12bfc069 | 63 | fHistMass(0), |
a273622f | 64 | fHistNEvents(0), |
12bfc069 | 65 | fVHF(0) |
3a219f60 | 66 | { |
9d6d35b0 | 67 | // Standard constructor |
3a219f60 | 68 | |
69 | // Output slot #1 writes into a TList container | |
70 | DefineOutput(1,TList::Class()); //My private output | |
10bdd1ae | 71 | // Output slot #2 writes into a TNtuple container |
72 | DefineOutput(2,TNtuple::Class()); //My private output | |
3a219f60 | 73 | } |
74 | ||
75 | //________________________________________________________________________ | |
76 | AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF() | |
77 | { | |
78 | // Destructor | |
79 | if (fOutput) { | |
80 | delete fOutput; | |
81 | fOutput = 0; | |
82 | } | |
12bfc069 | 83 | if (fVHF) { |
84 | delete fVHF; | |
85 | fVHF = 0; | |
86 | } | |
3a219f60 | 87 | } |
88 | ||
89 | //________________________________________________________________________ | |
90 | void AliAnalysisTaskSECompareHF::Init() | |
91 | { | |
92 | // Initialization | |
93 | ||
94 | if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n"); | |
9d6d35b0 | 95 | |
12bfc069 | 96 | gROOT->LoadMacro("ConfigVertexingHF.C"); |
97 | ||
98 | fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()"); | |
99 | fVHF->PrintStatus(); | |
9d6d35b0 | 100 | |
3a219f60 | 101 | return; |
102 | } | |
103 | ||
104 | //________________________________________________________________________ | |
105 | void AliAnalysisTaskSECompareHF::UserCreateOutputObjects() | |
106 | { | |
107 | // Create the output container | |
108 | // | |
109 | if(fDebug > 1) printf("AnalysisTaskSECompareHF::UserCreateOutputObjects() \n"); | |
110 | ||
111 | // Several histograms are more conveniently managed in a TList | |
112 | fOutput = new TList(); | |
113 | fOutput->SetOwner(); | |
114 | ||
115 | fHistMass = new TH1F("fHistMass", "D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965); | |
116 | fHistMass->Sumw2(); | |
117 | fHistMass->SetMinimum(0); | |
118 | fOutput->Add(fHistMass); | |
119 | ||
a273622f | 120 | fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5); |
121 | fHistNEvents->Sumw2(); | |
122 | fHistNEvents->SetMinimum(0); | |
123 | fOutput->Add(fHistNEvents); | |
124 | ||
2a5c77c6 | 125 | OpenFile(2); // 2 is the slot number of the ntuple |
607719af | 126 | fNtupleCmp = new TNtuple("fNtupleCmp","Charm comparison","pdg:nprongs:VxRec:VxTrue:ErrVx:VyRec:VyTrue:ErrVy:VzRec:VzTrue:ErrVz:Chi2toNDF:PtRec:Mrec:CPta:Prodd0"); |
3a219f60 | 127 | |
128 | return; | |
129 | } | |
130 | ||
131 | //________________________________________________________________________ | |
132 | void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/) | |
133 | { | |
134 | // Execute analysis for current event: | |
135 | // heavy flavor candidates association to MC truth | |
9d6d35b0 | 136 | |
3a219f60 | 137 | |
138 | AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent()); | |
9d6d35b0 | 139 | |
a273622f | 140 | fHistNEvents->Fill(0); // count event |
141 | // Post the data already here | |
142 | PostData(1,fOutput); | |
2a041947 | 143 | |
144 | // load HF vertices | |
145 | TClonesArray *inputArrayVertices = | |
146 | (TClonesArray*)aod->GetList()->FindObject("VerticesHF"); | |
147 | if(!inputArrayVertices) { | |
148 | printf("AliAnalysisTaskSECompareHF::UserExec: Vertices branch not found!\n"); | |
149 | return; | |
150 | } | |
151 | ||
3a219f60 | 152 | // load D0->Kpi candidates |
153 | TClonesArray *inputArrayD0toKpi = | |
154 | (TClonesArray*)aod->GetList()->FindObject("D0toKpi"); | |
155 | if(!inputArrayD0toKpi) { | |
156 | printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n"); | |
157 | return; | |
158 | } | |
159 | ||
9d6d35b0 | 160 | |
b3999999 | 161 | // load D*+ candidates |
162 | TClonesArray *inputArrayDstar = | |
163 | (TClonesArray*)aod->GetList()->FindObject("Dstar"); | |
164 | if(!inputArrayDstar) { | |
165 | printf("AliAnalysisTaskSECompareHF::UserExec: Dstar branch not found!\n"); | |
166 | return; | |
167 | } | |
9d6d35b0 | 168 | |
b3999999 | 169 | |
3a219f60 | 170 | // AOD primary vertex |
171 | AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex(); | |
172 | //vtx1->Print(); | |
173 | ||
174 | // load MC particles | |
175 | TClonesArray *mcArray = | |
176 | (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()); | |
177 | if(!mcArray) { | |
178 | printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n"); | |
179 | return; | |
180 | } | |
181 | ||
182 | // load MC header | |
183 | AliAODMCHeader *mcHeader = | |
184 | (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()); | |
185 | if(!mcHeader) { | |
186 | printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n"); | |
187 | return; | |
188 | } | |
189 | ||
2a041947 | 190 | Int_t nprongs,lab,okD0,okD0bar,pdg; |
191 | Bool_t unsetvtx; | |
90426506 | 192 | Double_t invmass,posRec[3],posTrue[3],covRec[6],errx,erry,errz; |
952f1e33 | 193 | AliAODRecoDecayHF2Prong *d2=0; |
194 | AliAODRecoDecayHF3Prong *d3=0; | |
195 | AliAODRecoDecayHF4Prong *d4=0; | |
2a041947 | 196 | |
197 | // loop over vertices | |
198 | Int_t nVertices = inputArrayVertices->GetEntriesFast(); | |
3b1598da | 199 | if(fDebug>1) printf("Number of vertices: %d\n",nVertices); |
2a041947 | 200 | |
201 | for (Int_t iVtx = 0; iVtx < nVertices; iVtx++) { | |
202 | AliAODVertex *vtx = (AliAODVertex*)inputArrayVertices->UncheckedAt(iVtx); | |
203 | ||
90426506 | 204 | vtx->GetXYZ(posRec); |
205 | vtx->GetCovarianceMatrix(covRec); | |
206 | errx=1.; erry=1.; errz=1.; | |
207 | if(covRec[0]>0) errx = TMath::Sqrt(covRec[0]); | |
208 | if(covRec[2]>0) erry = TMath::Sqrt(covRec[2]); | |
209 | if(covRec[5]>0) errz = TMath::Sqrt(covRec[5]); | |
210 | ||
211 | ||
2a041947 | 212 | // get parent AliAODRecoDecayHF |
213 | nprongs= vtx->GetNDaughters(); | |
607719af | 214 | // check that the daughters have kITSrefit and kTPCrefit |
215 | Bool_t allDgOK=kTRUE; | |
216 | for(Int_t idg=0; idg<nprongs; idg++) { | |
217 | AliAODTrack *track = (AliAODTrack*)vtx->GetDaughter(idg); | |
218 | if(!(track->GetStatus()&AliESDtrack::kITSrefit)) allDgOK=kFALSE; | |
219 | if(!(track->GetStatus()&AliESDtrack::kTPCrefit)) allDgOK=kFALSE; | |
220 | } | |
221 | if(!allDgOK) continue; | |
222 | ||
2a041947 | 223 | |
224 | switch(nprongs) { | |
225 | case 2: // look for D0->Kpi | |
952f1e33 | 226 | d2 = (AliAODRecoDecayHF2Prong*)vtx->GetParent(); |
27917274 | 227 | if(d2->IsLikeSign()) continue; |
a273622f | 228 | if(d2->Charge() != 0) continue; // these are D* |
2a041947 | 229 | lab = d2->MatchToMC(421,mcArray); |
230 | if(lab>=0) { | |
231 | unsetvtx=kFALSE; | |
232 | if(!d2->GetOwnPrimaryVtx()) { | |
233 | d2->SetOwnPrimaryVtx(vtx1); // needed to compute all variables | |
234 | unsetvtx=kTRUE; | |
235 | } | |
236 | okD0=0; okD0bar=0; | |
607719af | 237 | if(d2->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) { // beware: cuts may bias the resolution! |
238 | AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab); | |
239 | pdg = dMC->GetPdgCode(); | |
240 | invmass = (pdg==421 ? d2->InvMassD0() : d2->InvMassD0bar()); | |
607719af | 241 | // get a daughter for true pos of decay vertex |
242 | AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0)); | |
caf4ca8b | 243 | AliAODMCParticle *dg1MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(1)); |
607719af | 244 | dg0MC->XvYvZv(posTrue); |
caf4ca8b | 245 | // check that the D0 decayed in two prong (this is not necessary, because the check is already done, implicitely, in MatchToMC) |
246 | if(TMath::Abs(dMC->GetDaughter(0)-dMC->GetDaughter(1))!=1) continue; | |
247 | // check that the D0 decayed in K-pi+ | |
248 | if(!(TMath::Abs(dg0MC->GetPdgCode())==321 && TMath::Abs(dg1MC->GetPdgCode())==211) && | |
249 | !(TMath::Abs(dg0MC->GetPdgCode())==211 && TMath::Abs(dg1MC->GetPdgCode())==321)) continue; | |
250 | fHistMass->Fill(invmass); | |
251 | // Post the data already here | |
252 | PostData(1,fOutput); | |
607719af | 253 | Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs, |
254 | (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx, | |
255 | (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry, | |
256 | (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz, | |
257 | (Float_t)d2->GetReducedChi2(),(Float_t)d2->Pt(),(Float_t)invmass, | |
258 | (Float_t)d2->CosPointingAngle(),(Float_t)d2->Prodd0d0()}; | |
259 | fNtupleCmp->Fill(tmp); | |
260 | PostData(2,fNtupleCmp); | |
261 | } | |
2a041947 | 262 | if(unsetvtx) d2->UnsetOwnPrimaryVtx(); |
3a219f60 | 263 | } |
2a041947 | 264 | break; |
265 | case 3: // look for D+ | |
952f1e33 | 266 | d3 = (AliAODRecoDecayHF3Prong*)vtx->GetParent(); |
27917274 | 267 | if(d3->IsLikeSign()) continue; |
2a041947 | 268 | lab = d3->MatchToMC(411,mcArray); |
269 | if(lab>=0) { | |
270 | unsetvtx=kFALSE; | |
271 | if(!d3->GetOwnPrimaryVtx()) { | |
272 | d3->SetOwnPrimaryVtx(vtx1); // needed to compute all variables | |
273 | unsetvtx=kTRUE; | |
274 | } | |
90426506 | 275 | //if(d3->SelectDplus(fVHF->GetDplusCuts())) { |
276 | AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab); | |
277 | pdg = dMC->GetPdgCode(); | |
278 | invmass = d3->InvMassDplus(); | |
279 | // get a daughter for true pos of decay vertex | |
280 | AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0)); | |
281 | dg0MC->XvYvZv(posTrue); | |
caf4ca8b | 282 | // check that the D+ decayed in K-pi+pi+ |
283 | AliAODTrack *dg0 = (AliAODTrack*)d3->GetDaughter(0); | |
284 | AliAODTrack *dg1 = (AliAODTrack*)d3->GetDaughter(1); | |
285 | AliAODTrack *dg2 = (AliAODTrack*)d3->GetDaughter(2); | |
286 | dg0MC = (AliAODMCParticle*)mcArray->At(TMath::Abs(dg0->GetLabel())); | |
287 | AliAODMCParticle *dg1MC = (AliAODMCParticle*)mcArray->At(TMath::Abs(dg1->GetLabel())); | |
288 | AliAODMCParticle *dg2MC = (AliAODMCParticle*)mcArray->At(TMath::Abs(dg2->GetLabel())); | |
289 | ||
290 | if(!(TMath::Abs(dg0MC->GetPdgCode())==321 && TMath::Abs(dg1MC->GetPdgCode())==211 && TMath::Abs(dg2MC->GetPdgCode())==211) && | |
291 | !(TMath::Abs(dg0MC->GetPdgCode())==211 && TMath::Abs(dg1MC->GetPdgCode())==321 && TMath::Abs(dg2MC->GetPdgCode())==211) && | |
292 | !(TMath::Abs(dg0MC->GetPdgCode())==211 && TMath::Abs(dg1MC->GetPdgCode())==211 && TMath::Abs(dg2MC->GetPdgCode())==321)) continue; | |
607719af | 293 | Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs, |
294 | (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx, | |
295 | (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry, | |
296 | (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz, | |
297 | (Float_t)d3->GetReducedChi2(),(Float_t)d3->Pt(),(Float_t)invmass, | |
298 | (Float_t)d3->CosPointingAngle(),(Float_t)(d3->Getd0Prong(0)*d3->Getd0Prong(1)*d3->Getd0Prong(2))}; | |
299 | fNtupleCmp->Fill(tmp); | |
a7c1934a | 300 | PostData(2,fNtupleCmp); |
90426506 | 301 | //} |
2a041947 | 302 | if(unsetvtx) d3->UnsetOwnPrimaryVtx(); |
303 | } | |
304 | break; | |
305 | case 4: // look for D0->Kpipipi | |
952f1e33 | 306 | d4 = (AliAODRecoDecayHF4Prong*)vtx->GetParent(); |
27917274 | 307 | if(d4->IsLikeSign()) continue; |
2a041947 | 308 | lab = d4->MatchToMC(421,mcArray); |
309 | if(lab>=0) { | |
310 | unsetvtx=kFALSE; | |
311 | if(!d4->GetOwnPrimaryVtx()) { | |
312 | d4->SetOwnPrimaryVtx(vtx1); // needed to compute all variables | |
313 | unsetvtx=kTRUE; | |
314 | } | |
315 | okD0=0; okD0bar=0; | |
90426506 | 316 | //if(d4->SelectD0(fVHF->GetD0to4ProngsCuts(),okD0,okD0bar)) { |
317 | AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab); | |
318 | pdg = dMC->GetPdgCode(); | |
319 | //invmass = (pdg==421 ? d->InvMassD0() : d->InvMassD0bar()); | |
320 | invmass = 10.; // dummy | |
321 | // get a daughter for true pos of decay vertex | |
322 | AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0)); | |
323 | dg0MC->XvYvZv(posTrue); | |
607719af | 324 | Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs, |
325 | (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx, | |
326 | (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry, | |
327 | (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz, | |
328 | (Float_t)d4->GetReducedChi2(),(Float_t)d4->Pt(),(Float_t)invmass, | |
329 | (Float_t)d4->CosPointingAngle(),(Float_t)(d4->Getd0Prong(0)*d4->Getd0Prong(1)*d4->Getd0Prong(2)*d4->Getd0Prong(3))}; | |
330 | fNtupleCmp->Fill(tmp); | |
a7c1934a | 331 | PostData(2,fNtupleCmp); |
90426506 | 332 | //} |
2a041947 | 333 | if(unsetvtx) d4->UnsetOwnPrimaryVtx(); |
334 | } | |
335 | break; | |
3a219f60 | 336 | } |
3a219f60 | 337 | |
2a041947 | 338 | } // end loop on vertices |
3a219f60 | 339 | |
9d6d35b0 | 340 | |
b3999999 | 341 | // loop over D*+ candidates |
4dd17b11 | 342 | /* |
b3999999 | 343 | for (Int_t iDstar = 0; iDstar < inputArrayDstar->GetEntries(); iDstar++) { |
344 | AliAODRecoCascadeHF *c = (AliAODRecoCascadeHF*)inputArrayDstar->UncheckedAt(iDstar); | |
345 | Int_t labDstar = c->MatchToMC(413,421,mcArray); | |
346 | if(labDstar>=0) printf("GOOD MATCH FOR D*+\n"); | |
347 | } | |
4dd17b11 | 348 | */ |
2a041947 | 349 | |
4dd17b11 | 350 | // Post the data already here |
351 | PostData(1,fOutput); | |
2a041947 | 352 | PostData(2,fNtupleCmp); |
b3999999 | 353 | |
3a219f60 | 354 | return; |
355 | } | |
3a219f60 | 356 | //________________________________________________________________________ |
357 | void AliAnalysisTaskSECompareHF::Terminate(Option_t */*option*/) | |
358 | { | |
359 | // Terminate analysis | |
360 | // | |
361 | if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n"); | |
362 | ||
363 | fOutput = dynamic_cast<TList*> (GetOutputData(1)); | |
364 | if (!fOutput) { | |
365 | printf("ERROR: fOutput not available\n"); | |
366 | return; | |
367 | } | |
368 | ||
3a219f60 | 369 | fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass")); |
a273622f | 370 | fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents")); |
3a219f60 | 371 | |
27917274 | 372 | //fNtupleCmp = dynamic_cast<TNtuple*> (GetOutputData(2)); |
10bdd1ae | 373 | |
3a219f60 | 374 | return; |
375 | } | |
376 |