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