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