A little task for checking the c*tau of the strange particles
[u/mrichter/AliRoot.git] / PWG1 / global / AliAnalysisTaskCTau.cxx
1 #include <TCanvas.h>
2 #include <TTree.h>
3 #include <TFile.h>
4 #include <TH1F.h>
5 #include <TPDGCode.h>
6 #include <TDatabasePDG.h>
7 #include <TParticlePDG.h>
8 #include <TParticle.h>
9 #include <TROOT.h>
10
11 #include "AliESDEvent.h"
12 #include "AliESDv0.h"
13 #include "AliESDcascade.h"
14
15 #include "AliMCEvent.h"
16 #include "AliStack.h"
17
18 #include "AliAnalysisTaskCTau.h"
19
20 extern TROOT *gROOT;
21
22 //
23 //  This is a little task for checking the c*tau of the strange particles 
24 //
25
26 AliAnalysisTaskCTau::AliAnalysisTaskCTau(const char *name) :
27 AliAnalysisTaskSE(name),
28 fOutput(0),
29 fESD(0), 
30 fK0s(0),
31 fK0sMC(0),
32 fLambdas(0),
33 fLambdasMC(0),
34 fLambdaBars(0),
35 fLambdaBarsMC(0),
36 fXis(0),
37 fXisMC(0),
38 fMass(0),
39 fMassMC(0)
40 {
41   // Constructor. Initialization of pointers
42   DefineOutput(1, TList::Class());
43 }
44
45 void AliAnalysisTaskCTau::UserCreateOutputObjects()
46 {
47   fOutput = new TList(); 
48   fOutput->SetOwner();
49
50   fK0s     = new TH1F("K0s", "c\\tau for K^{0}_{s}", 100, 0., 25.);//10*2.684);
51   fK0s->GetXaxis()->SetTitle("L\\frac{m}{p} [cm]"); 
52   fOutput->Add(fK0s);
53   fK0sMC   = new TH1F("K0sMC", "c\\tau for K^{0}_{s}, MC", 100, 0., 25.);
54   fK0sMC->GetXaxis()->SetTitle("L\\frac{m}{p} [cm]"); 
55   fOutput->Add(fK0sMC);
56
57   fLambdas = new TH1F("Lambdas","c\\tau for \\Lambda", 100, 0.,80.);//10*7.89);
58   fLambdas->GetXaxis()->SetTitle("L\\frac{m}{p} [cm]"); 
59    //fLambdas->Sumw2();
60   fOutput->Add(fLambdas);
61   fLambdasMC = new TH1F("LambdasMC","c\\tau for \\Lambda, MC", 100, 0.,80.);
62   fLambdasMC->GetXaxis()->SetTitle("L\\frac{m}{p} [cm]"); 
63    //fLambdasMC->Sumw2();
64   fOutput->Add(fLambdasMC);
65
66   fLambdaBars = new TH1F("LambdaBars","c\\tau for \\bar{\\Lambda}",100,0.,80.);
67   fLambdaBars->GetXaxis()->SetTitle("L\\frac{m}{p} [cm]"); 
68    //fLambdaBars->Sumw2();
69   fOutput->Add(fLambdaBars);
70   fLambdaBarsMC=new TH1F("LambdaBarsMC","c\\tau for \\bar{\\Lambda}, MC",100,0.,80.);
71   fLambdaBarsMC->GetXaxis()->SetTitle("L\\frac{m}{p} [cm]"); 
72    //fLambdaBarsMC->Sumw2();
73   fOutput->Add(fLambdaBarsMC);
74
75   fXis = new TH1F("Xis","c\\tau for \\Xi^{-} + \\Xi^{+}",50,0.,30.);//10*4.91)
76   fXis->GetXaxis()->SetTitle("L\\frac{m}{p} [cm]"); 
77   fOutput->Add(fXis);
78   fXisMC = new TH1F("XisMC","c\\tau for \\Xi^{-} + \\Xi^{+}, MC",50,0.,30.);
79   fXisMC->GetXaxis()->SetTitle("L\\frac{m}{p} [cm]"); 
80   fOutput->Add(fXisMC);
81
82    //fMass = new TH1F("Mass", "Mass", 50, 0.4477, 0.5477);
83   fMass = new TH1F("Mass", "Mass", 50, 1.07, 1.16);
84    //fMass = new TH1F("Mass", "Mass", 50, 1.271, 1.371);
85   fOutput->Add(fMass);
86   fMassMC = new TH1F("MassMC", "Mass, MC", 50, 1.07, 1.16);
87   fOutput->Add(fMassMC);
88
89 }
90
91 void AliAnalysisTaskCTau::UserExec(Option_t *)
92 {
93
94   fESD=(AliESDEvent *)InputEvent();
95
96   if (!fESD) {
97     Printf("ERROR: fESD not available");
98     return;
99   }
100
101   //Char_t *tname="CINT1B-ABCE-NOPF-ALL";
102   //Int_t ok=fESD->IsTriggerClassFired(tname);
103   //if (!ok) return;
104   
105   const AliESDVertex *vtx=fESD->GetPrimaryVertexSPD();
106   //if (!vtx->GetStatus()) {
107   //   vtx=fESD->GetPrimaryVertexTracks();
108   //   if (!vtx->GetStatus()) return;
109   //}
110   Double_t xv=vtx->GetXv(), yv=vtx->GetYv(), zv=vtx->GetZv();
111
112
113   //+++++++ MC
114   AliStack *stack = 0x0;
115   AliMCEvent *mcEvent = MCEvent();
116   if (mcEvent) {
117      stack = mcEvent->Stack();
118      if (!stack) {
119         Printf("ERROR: stack not available");
120         return;
121      }
122
123      const AliVVertex *mcVtx=mcEvent->GetPrimaryVertex();
124      Double_t mcXv=mcVtx->GetX(), mcYv=mcVtx->GetY(), mcZv=mcVtx->GetZ();
125      Int_t ntrk=stack->GetNtrack();
126      while (ntrk--) {
127        TParticle *p0=stack->Particle(ntrk);
128        Int_t code=p0->GetPdgCode();
129        if (code != kK0Short)
130        if (code != kLambda0)
131        if (code != kLambda0Bar) continue;
132
133        Int_t plab=p0->GetFirstDaughter(), nlab=p0->GetLastDaughter();
134        if (nlab==plab) continue;
135        if (nlab<0) continue;
136        if (plab<0) continue;
137        TParticle *part = stack->Particle(plab);
138        if (!part) continue;
139        TParticlePDG *partPDG = part->GetPDG();
140        if (!partPDG) continue;
141        Double_t charge=partPDG->Charge();
142        if (charge == 0.) continue;
143        if (charge < 0.) {
144           Int_t i=plab; plab=nlab; nlab=i;
145        }
146   
147        Double_t pz=p0->Pz(), pt=p0->Pt(), p=p0->P();
148        if (pt<0.1) continue;
149        if (TMath::Abs(pz/pt)>0.7) continue;
150        if (p<0.9) continue;
151     
152        Double_t x=p0->Vx(), y=p0->Vy(), z=p0->Vz();
153        Double_t dx=mcXv-x, dy=mcYv-y, dz=mcZv-z;
154        Double_t l=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
155
156        if (l > 0.1) continue; // secondary V0
157
158        x=part->Vx(); y=part->Vy(); z=part->Vz();
159        dx=mcXv-x; dy=mcYv-y; dz=mcZv-z;
160        l=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
161
162        if (code == kK0Short) {
163           Double_t m=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
164           Double_t ct=l*m/p;    
165           fK0sMC->Fill(ct);
166        }
167        if (code == kLambda0) {
168           Double_t m=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
169           Double_t ct=l*m/p;    
170           fLambdasMC->Fill(ct);
171        }
172        if (code == kLambda0Bar) {
173          Double_t m=TDatabasePDG::Instance()->GetParticle(kLambda0Bar)->Mass();
174          Double_t ct=l*m/p;    
175          fLambdaBarsMC->Fill(ct);
176        }
177      }
178   }
179
180   Int_t nv0 = fESD->GetNumberOfV0s();
181   while (nv0--) {
182       AliESDv0 *v0=fESD->GetV0(nv0);
183
184       if (v0->GetOnFlyStatus()) continue;
185
186       Int_t nidx=TMath::Abs(v0->GetNindex());
187       AliESDtrack *ntrack=fESD->GetTrack(nidx);
188       if (!ntrack->IsOn(AliESDtrack::kTPCrefit)) continue;
189       //if (!ntrack->HasPointOnITSLayer(5)) continue;
190
191       Int_t pidx=TMath::Abs(v0->GetPindex());
192       AliESDtrack *ptrack=fESD->GetTrack(pidx);
193       if (!ptrack->IsOn(AliESDtrack::kTPCrefit)) continue;
194       //if (!ptrack->HasPointOnITSLayer(5)) continue;
195
196
197       Double_t pz=v0->Pz(),pt=v0->Pt(),p=v0->P();
198       if (pt<0.1) continue;
199       if (TMath::Abs(pz/pt)>0.7) continue;
200       if (p<0.9) continue;
201
202       Double_t cpa=v0->GetV0CosineOfPointingAngle();
203       if (cpa<0.995) continue;
204       
205
206       Double_t x,y,z; v0->GetXYZ(x,y,z);
207       Double_t dx=x-xv, dy=y-yv, dz=z-zv;
208       //Double_t r=TMath::Sqrt(x*x + y*y);
209       //if (r > 3.) continue;
210       Double_t l=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
211
212       v0->ChangeMassHypothesis(kK0Short);
213       Double_t mass=v0->GetEffMass();
214       Double_t m=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
215       if (TMath::Abs(m-mass)<0.016) {
216          Double_t ct=l*m/p;
217          fK0s->Fill(ct);
218       }
219       
220       //if (ntrack->GetP()<0.7) continue;
221
222       v0->ChangeMassHypothesis(kLambda0);
223       mass=v0->GetEffMass();
224       m=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
225       if (TMath::Abs(m-mass)<0.006) {
226          Double_t ct=l*m/p;
227          fLambdas->Fill(ct);
228       }
229
230       v0->ChangeMassHypothesis(kLambda0Bar);
231       mass=v0->GetEffMass();
232       m=TDatabasePDG::Instance()->GetParticle(kLambda0Bar)->Mass();
233       if (TMath::Abs(m-mass)<0.006) {
234          Double_t ct=l*m/p;
235          fLambdaBars->Fill(ct);
236       }
237       fMass->Fill(mass);
238
239
240       //+++++++ MC
241       Int_t nlab=TMath::Abs(ntrack->GetLabel());
242       Int_t plab=TMath::Abs(ptrack->GetLabel());
243       
244       if ((plab-nlab) != 1)
245         if ((nlab-plab) != 1) continue;
246
247       if (!stack) continue;
248       TParticle *np=stack->Particle(nlab);
249       Int_t i0=np->GetFirstMother();
250       TParticle *p0=stack->Particle(i0);
251       if (!p0) continue;
252       if (p0->GetPdgCode() == kLambda0Bar) fMassMC->Fill(mass);
253
254   }
255
256   Double_t kine0;
257   Int_t ncs=fESD->GetNumberOfCascades();
258   for (Int_t i=0; i<ncs; i++) {
259       AliESDcascade *cs=fESD->GetCascade(i);
260
261       Int_t nidx=TMath::Abs(cs->GetNindex());
262       AliESDtrack *ntrack=fESD->GetTrack(nidx);
263       if (!ntrack->IsOn(AliESDtrack::kTPCrefit)) continue;
264
265       Int_t pidx=TMath::Abs(cs->GetPindex());
266       AliESDtrack *ptrack=fESD->GetTrack(pidx);
267       if (!ptrack->IsOn(AliESDtrack::kTPCrefit)) continue;
268
269       Int_t bidx=TMath::Abs(cs->GetBindex());
270       AliESDtrack *btrack=fESD->GetTrack(bidx);
271       if (!btrack->IsOn(AliESDtrack::kTPCrefit)) continue;
272
273       Int_t charge=cs->Charge();      
274       if (charge<0)      
275          cs->ChangeMassHypothesis(kine0,kXiMinus);
276       else
277          cs->ChangeMassHypothesis(kine0,kXiPlusBar);
278
279       Double_t mass=cs->GetEffMassXi();
280
281       if (mass>1.305)
282       if (mass<1.335) {
283         //TTree *tree=fChain->GetTree();
284         //TDirectory *dir=tree->GetDirectory();
285         //Int_t n=fESD->GetEventNumberInFile();
286         //Printf("%d %d Xi found !!! %s %d\n", ncs, charge, dir->GetName(), n);
287         Double_t p=cs->P();
288         Double_t m=TDatabasePDG::Instance()->GetParticle(kXiMinus)->Mass();
289         Double_t x,y,z; cs->GetXYZcascade(x,y,z);
290         Double_t dx=x-xv, dy=y-yv, dz=z-zv;
291         Double_t l=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
292
293         if (p<0.1) continue;
294
295         Double_t ct=l*m/p;
296         fXis->Fill(ct);
297       }
298
299   }
300   
301   PostData(1, fOutput);
302 }
303
304 void AliAnalysisTaskCTau::Terminate(Option_t *)
305 {
306    // The Terminate() function is the last function to be called during
307    // a query. It always runs on the client, it can be used to present
308    // the results graphically or save the results to file.
309   
310   fOutput=(TList*)GetOutputData(1);
311   if (!fOutput) {
312      Printf("ERROR: fOutput not available");
313      return;
314   }
315  
316   fMass = dynamic_cast<TH1F*>(fOutput->FindObject("Mass")) ; 
317   if (!fMass) {
318      Printf("ERROR: fMass not available");
319      return;
320   }
321   fLambdas = dynamic_cast<TH1F*>(fOutput->FindObject("Lambdas")) ; 
322   if (!fLambdas) {
323      Printf("ERROR: fLambdas not available");
324      return;
325   }
326   fLambdaBars = dynamic_cast<TH1F*>(fOutput->FindObject("LambdaBars")) ; 
327   if (!fLambdaBars) {
328      Printf("ERROR: fLambdaBars not available");
329      return;
330   }
331   fXis = dynamic_cast<TH1F*>(fOutput->FindObject("Xis")) ; 
332   if (!fXis) {
333      Printf("ERROR: fXis not available");
334      return;
335   }
336   
337   if (!gROOT->IsBatch()) {
338     TCanvas *c1 = new TCanvas("c1","Lambdas");
339     c1->SetFillColor(10);
340     c1->SetHighLightColor(10);
341     fLambdas->DrawCopy("E") ;
342     TCanvas *c2 = new TCanvas("c2","LambdaBar");
343     c2->SetFillColor(10);
344     c2->SetHighLightColor(10);
345     fLambdaBars->DrawCopy("E") ;
346   }
347 }