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