remove option C for Clear for trigger array for the moment, causes malloc error
[u/mrichter/AliRoot.git] / PWG2 / 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   fMinTPC(-1E6),
38   fMaxTPC( 1E6),
39   fOutList(0x0),
40   fHEvents(0x0),
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());
51   DefineOutput(2, TTree::Class());
52   DefineOutput(3, TList::Class());
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),
70   fOutList(0x0),
71   fHEvents(0x0),
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
113   if (fOutTree[0]) delete fOutTree[0];
114   if (fOutTree[1]) delete fOutTree[1];
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
134   fTOFmaker->LoadChannelMap((char *)"tofmap.root");
135   fTOFmaker->SetMaskOffChannel();
136   
137   // initialize random
138   gRandom->SetSeed(0);
139
140   // create output trees
141   OpenFile(1);
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);
162 }
163
164
165 void AliRsnAnalysisPhi900GeV::UserExec(Option_t *)
166 {
167 //
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
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
191     if (v->GetNContributors() < 1)
192     {
193       fHEvents->Fill(3);
194       PostData(3, fOutList);
195       return;
196     }
197   }
198
199   // if the Z position is larger than 10, skip this event
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);
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
216   ProcessESD(esd, v, tof[0], stack);
217   ProcessMC(stack);
218
219   PostData(3, fOutList);
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   }
328   
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
378       fOutTree[0]->Fill();
379     }
380   }
381
382   PostData(1, fOutTree[0]);
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
435       fOutTree[1]->Fill();
436     }
437   }
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   }
465 }
466
467
468 Double_t AliRsnAnalysisPhi900GeV::AlephBB(Double_t p, Double_t mass) const
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