]>
Commit | Line | Data |
---|---|---|
ca8d2905 | 1 | // |
2 | // Implementation file for implementation of data analysis aft 900 GeV | |
3 | // | |
4 | // Author: A. Pulvirenti | |
5 | // | |
6 | ||
d7417f01 | 7 | #include "Riostream.h" |
8 | ||
9 | #include "TH1.h" | |
b9bbd271 | 10 | #include "TTree.h" |
11 | #include "TParticle.h" | |
12 | #include "TRandom.h" | |
13 | #include "TLorentzVector.h" | |
14 | ||
15 | #include "AliLog.h" | |
16 | #include "AliESDEvent.h" | |
17 | #include "AliESDVertex.h" | |
18 | #include "AliESDtrack.h" | |
19 | #include "AliStack.h" | |
20 | #include "AliMCEvent.h" | |
21 | ||
22 | #include "AliRsnAnalysisPhi900GeV.h" | |
23 | ||
24 | ||
25 | AliRsnAnalysisPhi900GeV::AliRsnAnalysisPhi900GeV(const char *name) : | |
26 | AliAnalysisTaskSE(name), | |
27 | fUseMC(kFALSE), | |
28 | fPDG(0), | |
29 | fIM(0.0), | |
30 | fPt(0.0), | |
31 | fY(0.0), | |
32 | fEta(0.0), | |
33 | fDCAr(1E6), | |
34 | fDCAz(1E6), | |
35 | fChi2(1E6), | |
36 | fNTPC(0), | |
37 | fMinTPC(-1E6), | |
38 | fMaxTPC( 1E6), | |
d7417f01 | 39 | fOutList(0x0), |
40 | fHEvents(0x0), | |
b9bbd271 | 41 | fTOFESD(kFALSE), |
42 | fTOFSigma(210.0), | |
43 | fTOFmaker(0x0), | |
44 | fTOFSettings(AliRsnTOFT0maker::kNone) | |
45 | { | |
46 | // | |
47 | // Constructor | |
48 | // | |
49 | ||
50 | DefineOutput(1, TTree::Class()); | |
d7417f01 | 51 | DefineOutput(2, TTree::Class()); |
52 | DefineOutput(3, TList::Class()); | |
b9bbd271 | 53 | } |
54 | ||
55 | ||
56 | AliRsnAnalysisPhi900GeV::AliRsnAnalysisPhi900GeV(const AliRsnAnalysisPhi900GeV& copy) : | |
57 | AliAnalysisTaskSE(copy), | |
58 | fUseMC(copy.fUseMC), | |
59 | fPDG(0), | |
60 | fIM(0.0), | |
61 | fPt(0.0), | |
62 | fY(0.0), | |
63 | fEta(0.0), | |
64 | fDCAr(copy.fDCAr), | |
65 | fDCAz(copy.fDCAz), | |
66 | fChi2(copy.fChi2), | |
67 | fNTPC(copy.fNTPC), | |
68 | fMinTPC(copy.fMinTPC), | |
69 | fMaxTPC(copy.fMaxTPC), | |
d7417f01 | 70 | fOutList(0x0), |
71 | fHEvents(0x0), | |
b9bbd271 | 72 | fTOFESD(copy.fTOFESD), |
73 | fTOFSigma(copy.fTOFSigma), | |
74 | fTOFmaker(0x0), | |
75 | fTOFSettings(copy.fTOFSettings) | |
76 | { | |
77 | // | |
78 | // Copy constructor | |
79 | // | |
80 | } | |
81 | ||
82 | ||
83 | AliRsnAnalysisPhi900GeV& AliRsnAnalysisPhi900GeV::operator=(const AliRsnAnalysisPhi900GeV& copy) | |
84 | { | |
85 | // | |
86 | // Assignment operator | |
87 | // | |
88 | ||
89 | fUseMC = copy.fUseMC; | |
90 | ||
91 | fDCAr = copy.fDCAr; | |
92 | fDCAz = copy.fDCAz; | |
93 | fChi2 = copy.fChi2; | |
94 | fNTPC = copy.fNTPC; | |
95 | ||
96 | fMinTPC = copy.fMinTPC; | |
97 | fMaxTPC = copy.fMaxTPC; | |
98 | ||
99 | fTOFESD = copy.fTOFESD; | |
100 | fTOFSigma = copy.fTOFSigma; | |
101 | fTOFSettings = copy.fTOFSettings; | |
102 | ||
103 | return (*this); | |
104 | } | |
105 | ||
106 | ||
107 | AliRsnAnalysisPhi900GeV::~AliRsnAnalysisPhi900GeV() | |
108 | { | |
109 | // | |
110 | // Destructor | |
111 | // | |
112 | ||
d7417f01 | 113 | if (fOutTree[0]) delete fOutTree[0]; |
114 | if (fOutTree[1]) delete fOutTree[1]; | |
b9bbd271 | 115 | } |
116 | ||
117 | ||
118 | void AliRsnAnalysisPhi900GeV::UserCreateOutputObjects() | |
119 | { | |
120 | // | |
121 | // Create the output data container | |
122 | // | |
123 | ||
124 | // setup TOF maker | |
125 | fTOFmaker = new AliRsnTOFT0maker; | |
126 | fTOFmaker->SetTimeResolution(fTOFSigma * 1E-12); | |
127 | fTOFmaker->SetESDdata(fTOFESD); | |
128 | fTOFmaker->fSettings = fTOFSettings; | |
129 | AliInfo(Form("TOF sigma = %f", fTOFSigma)); | |
130 | AliInfo(Form("TOF ESD = %s", (fTOFESD ? "YES" : "NO"))); | |
131 | AliInfo(Form("TOF settings = %s", fTOFmaker->Settings().Data())); | |
132 | ||
133 | // load dead channel map | |
6256671b | 134 | fTOFmaker->LoadChannelMap((char *)"tofmap.root"); |
b9bbd271 | 135 | fTOFmaker->SetMaskOffChannel(); |
136 | ||
137 | // initialize random | |
138 | gRandom->SetSeed(0); | |
139 | ||
140 | // create output trees | |
141 | OpenFile(1); | |
d7417f01 | 142 | fOutTree[0] = new TTree("rsnTree", "Pairs"); |
143 | ||
144 | fOutTree[0]->Branch("pdg", &fPDG, "pdg/S"); | |
145 | fOutTree[0]->Branch("im" , &fIM , "im/F" ); | |
146 | fOutTree[0]->Branch("y" , &fY , "y/F" ); | |
147 | fOutTree[0]->Branch("pt" , &fPt , "pt/F" ); | |
148 | fOutTree[0]->Branch("eta", &fEta, "eta/F"); | |
149 | ||
150 | OpenFile(2); | |
151 | fOutTree[1] = new TTree("rsnTrue", "True pairs"); | |
152 | ||
153 | fOutTree[1]->Branch("im" , &fIM , "im/F" ); | |
154 | fOutTree[1]->Branch("y" , &fY , "y/F" ); | |
155 | fOutTree[1]->Branch("pt" , &fPt , "pt/F" ); | |
156 | fOutTree[1]->Branch("eta", &fEta, "eta/F"); | |
157 | ||
158 | OpenFile(3); | |
159 | fOutList = new TList; | |
160 | fHEvents = new TH1I("hEvents", "Event details", 4, 0, 4); | |
161 | fOutList->Add(fHEvents); | |
b9bbd271 | 162 | } |
163 | ||
164 | ||
165 | void AliRsnAnalysisPhi900GeV::UserExec(Option_t *) | |
166 | { | |
167 | // | |
d7417f01 | 168 | // Main execution function. |
169 | // Fills the fHEvents data member with the following legenda: | |
170 | // 0 -- event OK, prim vertex with tracks | |
171 | // 1 -- event OK, prim vertex with SPD | |
172 | // 2 -- event OK but vz large | |
173 | // 3 -- event bad | |
b9bbd271 | 174 | // |
175 | ||
176 | // retrieve ESD event and related stack (if available) | |
177 | AliESDEvent *esd = dynamic_cast<AliESDEvent*>(fInputEvent); | |
178 | AliStack *stack = (fMCEvent ? fMCEvent->Stack() : 0x0); | |
179 | ||
180 | // get the best primary vertex: | |
181 | // first try the one with tracks | |
182 | Int_t type = 0; | |
183 | const AliESDVertex *v = esd->GetPrimaryVertexTracks(); | |
184 | if(v->GetNContributors() < 1) | |
185 | { | |
186 | // if not good, try SPD vertex | |
187 | type = 1; | |
188 | v = esd->GetPrimaryVertexSPD(); | |
189 | ||
190 | // if this is not good skip this event | |
d7417f01 | 191 | if (v->GetNContributors() < 1) |
192 | { | |
193 | fHEvents->Fill(3); | |
194 | PostData(3, fOutList); | |
195 | return; | |
196 | } | |
b9bbd271 | 197 | } |
198 | ||
199 | // if the Z position is larger than 10, skip this event | |
d7417f01 | 200 | if (TMath::Abs(v->GetZv()) > 10.0) |
201 | { | |
202 | fHEvents->Fill(2); | |
203 | PostData(3, fOutList); | |
204 | return; | |
205 | } | |
206 | ||
207 | // use the type to fill the histogram | |
208 | fHEvents->Fill(type); | |
b9bbd271 | 209 | |
210 | // smear TOF times in case of MC | |
211 | if (stack) RemakeTOFtimeMC(esd); | |
212 | ||
213 | // get time zero for TOF | |
214 | Double_t *tof = fTOFmaker->RemakePID(esd); | |
215 | ||
d7417f01 | 216 | ProcessESD(esd, v, tof[0], stack); |
217 | ProcessMC(stack); | |
218 | ||
219 | PostData(3, fOutList); | |
b9bbd271 | 220 | } |
221 | ||
222 | ||
223 | void AliRsnAnalysisPhi900GeV::Terminate(Option_t *) | |
224 | { | |
225 | // | |
226 | // Terminate | |
227 | // | |
228 | } | |
229 | ||
230 | ||
231 | void AliRsnAnalysisPhi900GeV::ProcessESD | |
232 | (AliESDEvent *esd, const AliESDVertex *v, Double_t time0, AliStack *stack) | |
233 | { | |
234 | // | |
235 | // This function works with the ESD object | |
236 | // | |
237 | ||
238 | // prepare to look on all tracks to select the ones | |
239 | // which pass all the cuts | |
240 | Int_t ntracks = esd->GetNumberOfTracks(); | |
241 | TArrayI pos(ntracks); | |
242 | TArrayI neg(ntracks); | |
243 | ||
244 | // define fixed functions for TOF compatibility range | |
245 | Double_t a1 = 0.01, a2 = -0.03; | |
246 | Double_t b1 = 0.25, b2 = 0.25; | |
247 | Double_t c1 = 0.05, c2 = -0.03; | |
248 | Double_t ymax, ymin; | |
249 | ||
250 | // loop on all tracks | |
251 | Int_t i, charge, nSPD, npos = 0, nneg = 0; | |
252 | Float_t chi2, b[2], bCov[3]; | |
253 | Double_t tpc, bb, mom, tofTime, tofRef, tofRel, times[10]; | |
254 | Bool_t okTOF; | |
255 | for (i = 0; i < ntracks; i++) | |
256 | { | |
257 | AliESDtrack *track = esd->GetTrack(i); | |
258 | if (!track) continue; | |
259 | ||
260 | // skip if it has not the required flags | |
261 | if (!track->IsOn(AliESDtrack::kTPCin)) continue; | |
262 | if (!track->IsOn(AliESDtrack::kTPCrefit)) continue; | |
263 | if (!track->IsOn(AliESDtrack::kITSrefit)) continue; | |
264 | ||
265 | // skip if it has not the TPC inner wall projection | |
266 | if (!track->GetInnerParam()) continue; | |
267 | ||
268 | // skip kink daughters | |
269 | if ((Int_t)track->GetKinkIndex(0) > 0) continue; | |
270 | ||
271 | // check clusters in TPC | |
272 | if (track->GetTPCclusters(0) < fNTPC) continue; | |
273 | ||
274 | // check chi2 | |
275 | chi2 = (Float_t)track->GetTPCchi2(); | |
276 | chi2 /= (Float_t)track->GetTPCclusters(0); | |
277 | if (chi2 > fChi2) continue; | |
278 | ||
279 | // check that has at least 1 SPD cluster | |
280 | nSPD = 0; | |
281 | if (track->HasPointOnITSLayer(0)) nSPD++; | |
282 | if (track->HasPointOnITSLayer(1)) nSPD++; | |
283 | if (nSPD < 1) continue; | |
284 | ||
285 | // check primary by reverting to vertex | |
286 | // and checking DCA | |
287 | if (!track->RelateToVertex(v, esd->GetMagneticField(), kVeryBig)) continue; | |
288 | track->GetImpactParameters(b, bCov); | |
289 | if (b[0] > fDCAr) continue; | |
290 | if (b[1] > fDCAz) continue; | |
291 | ||
292 | // check TPC dE/dx | |
293 | AliExternalTrackParam trackIn(*track->GetInnerParam()); | |
294 | mom = trackIn.P(); | |
295 | tpc = (Double_t)track->GetTPCsignal(); | |
296 | bb = AlephBB(mom); | |
297 | tpc = (tpc - bb) / bb; | |
298 | if (tpc < fMinTPC || tpc > fMaxTPC) continue; | |
299 | ||
300 | // if possible, check TOF | |
301 | okTOF = kTRUE; | |
302 | if (track->IsOn(AliESDtrack::kTOFpid)) | |
303 | { | |
304 | mom = track->P(); | |
305 | if (mom <= 0.26) | |
306 | okTOF = kTRUE; | |
307 | else | |
308 | { | |
309 | track->GetIntegratedTimes(times); | |
310 | tofTime = (Double_t)track->GetTOFsignal() - time0; | |
311 | tofRef = times[AliPID::kKaon]; | |
312 | tofRel = (tofTime - tofRef) / tofRef; | |
313 | ymax = a1 / (mom - b1) + c1; | |
314 | ymin = a2 / (mom - b2) + c2; | |
315 | okTOF = (tofRel >= ymin && tofRel <= ymax); | |
316 | } | |
317 | } | |
318 | if (!okTOF) continue; | |
319 | ||
320 | // if we arrive here, all cuts were passed | |
321 | // and we add the track to one array depending on charge | |
322 | charge = (Int_t)track->Charge(); | |
323 | if (charge > 0) | |
324 | pos[npos++] = i; | |
325 | else if (charge < 0) | |
326 | neg[nneg++] = i; | |
327 | } | |
d7417f01 | 328 | |
b9bbd271 | 329 | // resize arrays accordingly |
330 | pos.Set(npos); | |
331 | neg.Set(nneg); | |
332 | ||
333 | // loop to compute invariant mass | |
334 | Int_t ip, in, lp, ln; | |
335 | AliPID pid; | |
336 | Double_t kmass = pid.ParticleMass(AliPID::kKaon); | |
337 | Double_t phimass = 1.019455; | |
338 | TParticle *partp = 0x0, *partn = 0x0; | |
339 | AliESDtrack *tp = 0x0, *tn = 0x0; | |
340 | TLorentzVector vp, vn, vsum, vref; | |
341 | for (ip = 0; ip < npos; ip++) | |
342 | { | |
343 | tp = esd->GetTrack(pos[ip]); | |
344 | lp = TMath::Abs(tp->GetLabel()); | |
345 | if (stack) partp = stack->Particle(lp); | |
346 | ||
347 | for (in = 0; in < nneg; in++) | |
348 | { | |
349 | if (pos[ip] == neg[in]) continue; | |
350 | tn = esd->GetTrack(neg[in]); | |
351 | ln = TMath::Abs(tn->GetLabel()); | |
352 | if (stack) partn = stack->Particle(ln); | |
353 | ||
354 | fPDG = 0; | |
355 | if (partp && partn) | |
356 | { | |
357 | if (partp->GetFirstMother() == partn->GetFirstMother()) | |
358 | { | |
359 | if (partp->GetFirstMother() > 0) | |
360 | { | |
361 | TParticle *mum = stack->Particle(partp->GetFirstMother()); | |
362 | fPDG = mum->GetPdgCode(); | |
363 | } | |
364 | } | |
365 | } | |
366 | fPDG = TMath::Abs(fPDG); | |
367 | ||
368 | vp.SetXYZM(tp->Px(), tp->Py(), tp->Pz(), kmass); | |
369 | vn.SetXYZM(tn->Px(), tn->Py(), tn->Pz(), kmass); | |
370 | vsum = vp + vn; | |
371 | vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass); | |
372 | ||
373 | fIM = (Float_t)vsum.M(); | |
374 | fPt = (Float_t)vsum.Perp(); | |
375 | fEta = (Float_t)vsum.Eta(); | |
376 | fY = (Float_t)vref.Rapidity(); | |
377 | ||
d7417f01 | 378 | fOutTree[0]->Fill(); |
b9bbd271 | 379 | } |
380 | } | |
d7417f01 | 381 | |
382 | PostData(1, fOutTree[0]); | |
b9bbd271 | 383 | } |
384 | ||
385 | ||
386 | void AliRsnAnalysisPhi900GeV::ProcessMC(AliStack *stack) | |
387 | { | |
388 | // | |
389 | // Function to process stack only | |
390 | // | |
391 | ||
392 | if (!stack) return; | |
393 | Int_t nPart = stack->GetNtrack(); | |
394 | ||
395 | // loop to compute invariant mass | |
396 | Int_t ip, in; | |
397 | AliPID pid; | |
398 | Double_t kmass = pid.ParticleMass(AliPID::kKaon); | |
399 | Double_t phimass = 1.019455; | |
400 | TParticle *partp = 0x0, *partn = 0x0; | |
401 | TLorentzVector vp, vn, vsum, vref; | |
402 | ||
403 | for (ip = 0; ip < nPart; ip++) | |
404 | { | |
405 | partp = stack->Particle(ip); | |
406 | if (partp->GetPdgCode() != 321) continue; | |
407 | ||
408 | for (in = 0; in < nPart; in++) | |
409 | { | |
410 | partn = stack->Particle(in); | |
411 | if (partn->GetPdgCode() != -321) continue; | |
412 | ||
413 | fPDG = 0; | |
414 | if (partp->GetFirstMother() == partn->GetFirstMother()) | |
415 | { | |
416 | if (partp->GetFirstMother() > 0) | |
417 | { | |
418 | TParticle *mum = stack->Particle(partp->GetFirstMother()); | |
419 | fPDG = mum->GetPdgCode(); | |
420 | } | |
421 | } | |
422 | fPDG = TMath::Abs(fPDG); | |
423 | if (fPDG != 333) continue; | |
424 | ||
425 | vp.SetXYZM(partp->Px(), partp->Py(), partp->Pz(), kmass); | |
426 | vn.SetXYZM(partn->Px(), partn->Py(), partn->Pz(), kmass); | |
427 | vsum = vp + vn; | |
428 | vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass); | |
429 | ||
430 | fIM = (Float_t)vsum.M(); | |
431 | fPt = (Float_t)vsum.Perp(); | |
432 | fEta = (Float_t)vsum.Eta(); | |
433 | fY = (Float_t)vref.Rapidity(); | |
434 | ||
d7417f01 | 435 | fOutTree[1]->Fill(); |
b9bbd271 | 436 | } |
437 | } | |
d7417f01 | 438 | |
439 | PostData(2, fOutTree[1]); | |
440 | } | |
441 | ||
442 | ||
443 | void AliRsnAnalysisPhi900GeV::SetTPCparams(Bool_t isMC) | |
444 | { | |
445 | // | |
446 | // Set TPC bethe-bloch parameters | |
447 | // | |
448 | ||
449 | if (!isMC) | |
450 | { | |
451 | fTPCpar[0] = 1.41543; | |
452 | fTPCpar[1] = 2.63394E1; | |
453 | fTPCpar[2] = 5.0411E-11; | |
454 | fTPCpar[3] = 2.12543; | |
455 | fTPCpar[4] = 4.88663; | |
456 | } | |
457 | else | |
458 | { | |
459 | fTPCpar[0] = 2.15898; | |
460 | fTPCpar[1] = 1.75295E1; | |
461 | fTPCpar[2] = 3.40030E-9; | |
462 | fTPCpar[3] = 1.96178; | |
463 | fTPCpar[4] = 3.91720; | |
464 | } | |
b9bbd271 | 465 | } |
466 | ||
467 | ||
ca8d2905 | 468 | Double_t AliRsnAnalysisPhi900GeV::AlephBB(Double_t p, Double_t mass) const |
b9bbd271 | 469 | { |
470 | // | |
471 | // Compute expected Bethe-Bloch for that momentum and mass | |
472 | // | |
473 | ||
474 | if (mass < 1E-6) return 0.0; | |
475 | ||
476 | Double_t aa, bb, out, beta, bg; | |
477 | ||
478 | bg = p / mass; | |
479 | beta = bg / TMath::Sqrt(1.0 + bg * bg); | |
480 | aa = TMath::Power(beta, fTPCpar[3]); | |
481 | bb = TMath::Power(1./bg, fTPCpar[4]); | |
482 | bb = TMath::Log(fTPCpar[2] + bb); | |
483 | out = (fTPCpar[1] - aa - bb) * fTPCpar[0]/aa; | |
484 | ||
485 | return out; | |
486 | } | |
487 | ||
488 | ||
489 | Double_t AliRsnAnalysisPhi900GeV::RemakeTOFtimeMC(AliESDEvent *& esd) | |
490 | { | |
491 | // | |
492 | // Smear initial time for TOF in order to reproduce data resolution | |
493 | // | |
494 | ||
495 | Double_t t0 = gRandom->Gaus(0,135.); //spread in ps | |
496 | Int_t ntracks = esd->GetNumberOfTracks(); | |
497 | Double_t sigmaStandard = 80.; // 80 ps from TOF | |
498 | ||
499 | while (ntracks--) | |
500 | { | |
501 | AliESDtrack *t = esd->GetTrack(ntracks); | |
502 | if ((t->GetStatus()&AliESDtrack::kTOFout) == 0 || (t->GetStatus()&AliESDtrack::kTIME)==0) continue; | |
503 | Double_t time = t->GetTOFsignal(); | |
504 | if(fTOFSigma > sigmaStandard) | |
505 | { | |
506 | Double_t sigmaAdded = TMath::Sqrt(fTOFSigma*fTOFSigma - sigmaStandard*sigmaStandard); | |
507 | Double_t timerandomtrack = gRandom->Gaus(0, sigmaAdded); //spread in ps | |
508 | time += timerandomtrack; | |
509 | } | |
510 | time += t0; | |
511 | t->SetTOFsignal(time); | |
512 | } | |
513 | return t0; | |
514 | } |