Fixed bug in the computation of dip angle value
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnAnalysisPhi900GeV.cxx
CommitLineData
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
25AliRsnAnalysisPhi900GeV::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
56AliRsnAnalysisPhi900GeV::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
83AliRsnAnalysisPhi900GeV& 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
107AliRsnAnalysisPhi900GeV::~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
118void 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
165void 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
223void AliRsnAnalysisPhi900GeV::Terminate(Option_t *)
224{
225//
226// Terminate
227//
228}
229
230
231void 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
386void 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
443void 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 468Double_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
489Double_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}