]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/cosmic/AliAnalysisTaskCosmic.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
make Kalman propagator more flexible for testing
[u/mrichter/AliRoot.git] / PWG1 / cosmic / AliAnalysisTaskCosmic.cxx
CommitLineData
c2174181 1#include "TChain.h"
2#include "TTree.h"
3#include "TH1F.h"
4#include "TH2F.h"
5#include "TProfile.h"
6#include "TList.h"
7#include "TTree.h"
8#include "TParticle.h"
9#include "TParticlePDG.h"
10
11
12#include "AliESDEvent.h"
13#include "AliESDtrack.h"
14#include "AliAnalysisTaskCosmic.h"
15#include "AliESDInputHandler.h"
16
17ClassImp(AliAnalysisTaskCosmic)
18
19//________________________________________________________________________
20AliAnalysisTaskCosmic::AliAnalysisTaskCosmic(const char *name)
21 : AliAnalysisTaskSE(name),
22 fHists(0),
23 fhDZvsZ(0),
24 fhDZvsPhi(0),
25 fhCh1Ch2(0),
26 fhPh1Ph2(0),
27 fhCl1Cl2G(0),
28 fhCl1Cl2B(0)
29{
30 //
31 // Constructor
32 //
33 for (Int_t i = 0; i < 6; i++) {
34 fhPt[i] = 0;
35 fhTheta[i] = 0;
36 fhPhi[i] = 0;
37 fhDPhi[i] = 0;
38 fhDTheta[i] = 0;
39 fhDZ[i] = 0;
40 fhDX[i] = 0;
41 fhDY[i] = 0;
42 fhDPt[i] = 0;
43 fhD1ovPt[i] = 0;
44 fpDPt[i] = 0;
45 fhDPtovPt[i] = 0;
46 fpDPtS[i] = 0;
47 }
48
49 DefineOutput(1, TList::Class());
50}
51
52
53//________________________________________________________________________
54void AliAnalysisTaskCosmic::UserCreateOutputObjects()
55{
56 // Create histograms
57 // Called once
58 TString ext[6] = {"PC", "NC", "PZ", "NZ", "_Good", "_Bad"};
59
60 char name[12];
61
62 fHists = new TList();
63
64 for (Int_t i = 0; i < 6; i++) {
65 // Pt
66 sprintf(name, "fhPt%2s", ext[i].Data());
67 fhPt[i] = new TH1F(name, " pT distribution", 800, 0., 200.);
68 fhPt[i]->SetXTitle("p_{T} [Gev]");
69 // Phi
70 sprintf(name, "fhPhi%2s", ext[i].Data());
71 fhPhi[i] = new TH1F(name, "Phi distribution", 62, 0., 2. * TMath::Pi());
72 fhPhi[i]->SetXTitle("#phi [rad]");
73 // Theta
74 sprintf(name, "fhTheta%2s", ext[i].Data());
75 fhTheta[i] = new TH1F(name, "Theta distribution", 62, 0., TMath::Pi());
76 fhTheta[i]->SetXTitle("#theta [rad]");
77 // Delta Phi
78 sprintf(name, "fhDPhi%2s", ext[i].Data());
79 fhDPhi[i] = new TH1F(name, "DeltaPhi distribution", 320, -0.4, 0.4);
80 fhDPhi[i]->SetXTitle("#Delta#phi [rad]");
81 // Delta Theta
82 sprintf(name, "fhDTheta%2s", ext[i].Data());
83 fhDTheta[i] = new TH1F(name, "DeltaTheta distribution", 320, -0.4, 0.4);
84 fhDTheta[i]->SetXTitle("#Delta#theta [rad]");
85 // Delta Z
86 sprintf(name, "fhDZ%2s", ext[i].Data());
87 fhDZ[i] = new TH1F(name, "DeltaZ distribution", 200, -10., 10.);
88 fhDZ[i]->SetXTitle("#DeltaZ [cm]");
89 // Delta X
90 sprintf(name, "fhDX%2s", ext[i].Data());
91 fhDX[i] = new TH1F(name, "DeltaX distribution", 200, -10., 10.);
92 fhDX[i]->SetXTitle("#DeltaX [cm]");
93 // Delta Y
94 sprintf(name, "fhDY%2s", ext[i].Data());
95 fhDY[i] = new TH1F(name, "DeltaY distribution", 200, -10, 10.);
96 fhDY[i]->SetXTitle("#DeltaY [cm]");
97 // Delta Pt
98 sprintf(name, "fhDPt%2s", ext[i].Data());
99 fhDPt[i] = new TH1F(name, "DeltaPt distribution", 200, -20., 20.);
100 fhDPt[i]->SetXTitle("#Delta p_{T} [GeV]");
101
102 // Delta 1/Pt
103 sprintf(name, "fhD1ovPt%2s", ext[i].Data());
104 fhD1ovPt[i] = new TH1F(name, "Delta 1/Pt distribution", 200, -1., 1.);
105 fhD1ovPt[i]->SetXTitle("#Delta 1/Pt");
106
107 // Delta Pt over Pt
108 sprintf(name, "fhDPtovPt%2s", ext[i].Data());
109 fhDPtovPt[i] = new TH1F(name, "DeltaPt/Pt distribution", 200, -2., 2.);
110 fhDPtovPt[i]->SetXTitle("#DeltaPt/Pt");
111
112
113
114 // Delta Pt/ Pt vs Pt
115 sprintf(name, "fpDPt%2s", ext[i].Data());
116 fpDPt[i] = new TProfile(name, "#Delta Pt / Pt", 20, 0., 20., -1, 1., "S");
117 fpDPt[i]->SetXTitle("p_{T} [GeV]");
118 fpDPt[i]->SetYTitle("#Delta 1/p_{T} [GeV^{-1}]");
119 // Delta Pt error
120 sprintf(name, "fpDPtS%2s", ext[i].Data());
121 fpDPtS[i] = new TProfile(name, "#Delta Pt / <sigma>", 20, 0., 20., 0., 10.);
122 fpDPtS[i]->SetXTitle("p_{T}");
123 fpDPtS[i]->SetYTitle("#Delta p_{T} / <#sigma_{p_{T}}>");
124
125
126 fHists->Add(fhPt[i]);
127 fHists->Add(fhPhi[i]);
128 fHists->Add(fhTheta[i]);
129 fHists->Add(fhDPhi[i]);
130 fHists->Add(fhDPt[i]);
131 fHists->Add(fhD1ovPt[i]);
132 fHists->Add(fhDPtovPt[i]);
133 fHists->Add(fhDTheta[i]);
134 fHists->Add(fhDZ[i]);
135 fHists->Add(fhDX[i]);
136 fHists->Add(fhDY[i]);
137 fHists->Add(fpDPt[i]);
138 fHists->Add(fpDPtS[i]);
139 }
140
141 fhDZvsZ = new TH2F("fhDZvsZ", "dz vs z", 500, -250., 250., 100, -10., 10.);
142 fhDZvsZ->SetXTitle("z_{in} * sign(z_{in}) * sign(z_{out}) [cm]");
143 fhDZvsZ->SetYTitle("#Delta z [cm]");
144
145
146 fhDZvsPhi = new TH2F("fhDZvsPhi", "dz vs phi", 36, 0., TMath::Pi(), 50, -2., 2.);
147
148 fhCh1Ch2 = new TH2F("fCh1Ch2", "ch1 vs ch2", 8, -2., 2., 8, -2., 2.);
149 fhPh1Ph2 = new TH2F("fPh1Ph2", "ph1 vs ph2", 128, 0., 2. * TMath::Pi(), 128, 0., 2. * TMath::Pi());
150 fhCl1Cl2G = new TH2F("fCl1Cl2G", "#cl vs #cl", 200, 0., 200., 200, 0., 200);
151 fhCl1Cl2B = new TH2F("fCl1Cl2B", "#cl vs #cl", 200, 0., 200., 200, 0., 200);
152
153 fHists->Add(fhDZvsZ);
154 fHists->Add(fhDZvsPhi);
155 fHists->Add(fhCh1Ch2);
156 fHists->Add(fhPh1Ph2);
157 fHists->Add(fhCl1Cl2G);
158 fHists->Add(fhCl1Cl2B);
159 fHists->SetOwner();
160
161}
162
163//________________________________________________________________________
164void AliAnalysisTaskCosmic::UserExec(Option_t *)
165{
166 // Main loop
167 // Called for each event
168
169 if (!fInputEvent) {
170 Printf("ERROR: fESD not available");
171 return;
172 }
173
174 AliESDEvent* esdE = (AliESDEvent*)fInputEvent;
86e75636 175
176 if (esdE->GetNumberOfTracks() != 2) {
177 PostData(1, fHists);
178 return;
179 }
180
c2174181 181
182 for (Int_t iTracks = 0; iTracks < esdE->GetNumberOfTracks(); iTracks++) {
183 AliESDtrack* track = esdE->GetTrack(iTracks);
184
185
186
187// const AliExternalTrackParam * track = trackG->GetTPCInnerParam();
188// if (!track) continue;
189 // Pt cut
190 Float_t pt = track->Pt();
191 Float_t charge = track->Charge();
192 Float_t phi = track->Phi();
193 if (phi > 0. && phi < TMath::Pi()) charge*=(-1.);
194
195
196
197 // TPC track
198 UInt_t status = track->GetStatus();
199 if (status&AliESDtrack::kTPCrefit ==0) continue;
200
201 Int_t nClustersTPC = track->GetTPCclusters(0);
202 if (nClustersTPC < 50) continue;
203
204 //
205
206 Float_t z = track->GetZ();
207 Float_t x = track->Xv();
208 Float_t y = track->Yv();
209 Float_t theta = track->Theta();
210
211
212 const AliExternalTrackParam * trackOut = track->GetOuterParam();
213 Float_t zOut = 0.;
214 if (trackOut)zOut = trackOut->Zv();
215
216
217 if (charge > 0) {
218 fhPt[kPosC] ->Fill(pt);
219 } else {
220 fhPt[kNegC] ->Fill(pt);
221 }
222
223 // if ((TMath::Abs(esdE->GetCurrentL3()) > 1.e-3) && (pt < 1. || pt > 50.)) continue;
224
225
226 //
227 if (charge > 0) {
228 fhPhi[kPosC] ->Fill(phi);
229 fhTheta[kPosC]->Fill(theta);
230 } else {
231 fhPhi[kNegC] ->Fill(phi);
232 fhTheta[kNegC]->Fill(theta);
233 }
234
235
236 if (z > 0) {
237 fhPt[kPosZ] ->Fill(pt);
238 fhPhi[kPosZ] ->Fill(phi);
239 fhTheta[kPosZ]->Fill(theta);
240 } else {
241 fhPt[kNegZ] ->Fill(pt);
242 fhPhi[kNegZ] ->Fill(phi);
243 fhTheta[kNegZ]->Fill(theta);
244 }
245
246
247
248 // Tracks coming from above
249 if (phi > 0. && phi < TMath::Pi()) {
250
251// printf("Track#1 %5d %5d %13.3f %13.3f \n", Int_t(Entry()), iTracks, phi, zOut);
252
253 Float_t dphiMin = 999.;
254 Float_t rMin = 999.;
255 Int_t jMin = -1;
256 // Search for a matching track (in dphi-dtheta)
257 for (Int_t jTracks = 0; jTracks < esdE->GetNumberOfTracks(); jTracks++) {
258 if (jTracks == iTracks) continue;
259
260 AliESDtrack* track2 = esdE->GetTrack(jTracks);
261
262 UInt_t status = track2->GetStatus();
263 if (status&AliESDtrack::kTPCrefit ==0) continue;
264
265 Int_t nClustersTPC2 = track2->GetTPCclusters(0);
266 if (nClustersTPC2 < 50) continue;
267 // if ((TMath::Abs(esdE->GetCurrentL3()) > 1.e-3) && (track2->Pt() < 1. || track2->Pt() > 50.)) continue;
268// const AliExternalTrackParam * track2 = trackG2->GetTPCInnerParam();
269// if (!track2) continue;
270
271 Float_t phi2 = track2->Phi() - TMath::Pi();
272 Float_t theta2 = TMath::Pi() - track2->Theta();
273 Float_t dphi = phi2 - phi;
274 Float_t dtheta = theta2 - theta;
275
276 if (dphi > TMath::Pi()) dphi -= 2. * TMath::Pi();
277 if (dphi < -TMath::Pi()) dphi += 2. * TMath::Pi();
278
279 Float_t dR = TMath::Sqrt(dphi * dphi + dtheta * dtheta);
280
281 if (dR < rMin) {
282 rMin = dR;
283 dphiMin = dphi;
284 jMin = jTracks;
285 }
286
287 } // tracks 2
288 if (jMin != -1) {
289 // we found a matching track candidate ...
290 AliESDtrack* track2 = esdE->GetTrack(jMin);
291 const AliExternalTrackParam * trackOut2 = track2->GetOuterParam();
292 Float_t zOut2 = 0.;
293 if (trackOut2) zOut2 = trackOut2->Zv();
294
295
296 Float_t theta2 = - track2->Theta() + TMath::Pi();
297 Float_t z2 = track2->GetZ();
298 Float_t x2 = track2->Xv();
299 Float_t y2 = track2->Yv();
300 Float_t charge2 = track2->Charge();
301 Float_t dz = z2 - z;
302 Float_t dx = x2 - x;
303 Float_t dy = y2 - y;
304 Float_t pt1 = track->Pt();
305 Float_t pt2 = track2->Pt();
306 Float_t dpt = pt2 - pt1;
307 Float_t d1pt = 1./pt2 - 1./pt1;
308
309 Float_t ptm = 0.5 * (pt1 + pt2);
310
311 Int_t nClustersTPC2 = track2->GetTPCclusters(0);
312
313 if (charge > 0.) {
314 fhDPhi[kPosC] ->Fill(dphiMin);
315 fhDTheta[kPosC]->Fill(theta2 - theta);
316 } else {
317 fhDPhi[kNegC] ->Fill(dphiMin);
318 fhDTheta[kNegC]->Fill(theta2 - theta);
319 }
320
321 if (z > 0.) {
322 fhDPhi[kPosZ] ->Fill(dphiMin);
323 fhDTheta[kPosZ]->Fill(theta2 - theta);
324 } else {
325 fhDPhi[kNegZ] ->Fill(dphiMin);
326 fhDTheta[kNegZ]->Fill(theta2 - theta);
327 }
328
329 if (TMath::Abs(dpt)/ptm > 0.5) {
330 fhDPhi[kBad] ->Fill(dphiMin);
331 fhDTheta[kBad]->Fill(theta2 - theta);
332 } else {
333 fhDPhi[kGood] ->Fill(dphiMin);
334 fhDTheta[kGood]->Fill(theta2 - theta);
335 }
336 // Good matches ...
337// if (TMath::Abs(rMin < 0.04) && (charge == charge2) && TMath::Abs(dz) < 0.5 && TMath::Abs(dy) && TMath::Abs(dx) < 0.5 )
338// if (TMath::Abs(rMin < 0.04) && (charge == charge2) && (zOut * zOut2) < 0.)
339 if (TMath::Abs(rMin < 0.04) && (charge == charge2))
340 {
341
342 if (TMath::Abs(dpt)/ptm < .1) {
343 fhCl1Cl2G->Fill(nClustersTPC, nClustersTPC2);
344 }
345
346 if (TMath::Abs(dpt)/ptm > .5) fhCl1Cl2B->Fill(nClustersTPC, nClustersTPC2);
347 fhPh1Ph2->Fill(track->Phi(), track2->Phi());
348
349
350 Double_t sigmaPt1 = TMath::Sqrt(track->GetSigma1Pt2());
351 Double_t sigmaPt2 = TMath::Sqrt(track->GetSigma1Pt2());
352 Double_t sigmaPt = 0.5 * TMath::Sqrt(sigmaPt1 * sigmaPt1 + sigmaPt2 * sigmaPt2);
353 if (TMath::Abs(dpt)/ptm > 0.2 && pt > 10. && charge > 0.)
354// printf("Track#2 %5d %5d %13.3f %13.3f %13.3f %13.3f\n", Entry(), jMin, track2->Phi(), dz, dpt, zOut2);
355 if (zOut != 0. && zOut2 != 0.) fhDZvsZ->Fill(TMath::Abs(zOut) * zOut/zOut2, dz);
356 if (zOut * zOut2 > 0. && zOut < 0) fhDZvsPhi->Fill(phi, dz);
357
358 fhCh1Ch2->Fill(charge, track2->Charge());
359
360
361 if (charge > 0.) {
362 fhDPt[kPosC]->Fill(dpt);
363 fhD1ovPt[kPosC]->Fill(d1pt);
364 fpDPt[kPosC]->Fill(ptm, d1pt);
365 if (ptm > 5. && ptm < 15.) fhDPtovPt[kPosC]->Fill(dpt/ptm);
366 fhDZ[kPosC]->Fill(dz);
367 fhDX[kPosC]->Fill(dx);
368 fhDY[kPosC]->Fill(dy);
369 fpDPtS[kPosC]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
370 } else {
371 fhDPt[kNegC]->Fill(dpt);
372 fhD1ovPt[kNegC]->Fill(d1pt);
373 fpDPt[kNegC]->Fill(ptm, d1pt);
374 if (ptm > 5. && ptm < 15.) fhDPtovPt[kNegC]->Fill(dpt/ptm);
375 fhDZ[kNegC]->Fill(dz);
376 fhDX[kNegC]->Fill(dx);
377 fhDY[kNegC]->Fill(dy);
378 fpDPtS[kNegC]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
379 }
380
381 if (z > 0.) {
382 fhDPt[kPosZ]->Fill(dpt);
383 fhD1ovPt[kPosZ]->Fill(d1pt);
384 fpDPt[kPosZ]->Fill(ptm, d1pt);
385 fhDZ[kPosZ]->Fill(dz);
386 fhDX[kPosZ]->Fill(dx);
387 fhDY[kPosZ]->Fill(dy);
388 fpDPtS[kPosZ]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
389 } else {
390 fhDPt[kNegZ]->Fill(dpt);
391 fhD1ovPt[kNegZ]->Fill(d1pt);
392 fpDPt[kNegZ]->Fill(ptm, d1pt);
393 fhDZ[kNegZ]->Fill(dz);
394 fhDX[kNegZ]->Fill(dx);
395 fhDY[kNegZ]->Fill(dy);
396 fpDPtS[kNegZ]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
397 }
398
399
400 if (TMath::Abs(dpt)/ptm > 0.5) {
401 fhDPt[kBad]->Fill(dpt);
402 fpDPt[kBad]->Fill(ptm, TMath::Abs(dpt)/ptm);
403 fhDZ[kBad]->Fill(dz);
404 fhDX[kBad]->Fill(dx);
405 fhDY[kBad]->Fill(dy);
406 fpDPtS[kBad]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
407 } else {
408 fhDPt[kGood]->Fill(dpt);
409 fpDPt[kGood]->Fill(ptm, TMath::Abs(dpt)/ptm);
410 fhDZ[kGood]->Fill(dz);
411 fhDX[kGood]->Fill(dx);
412 fhDY[kGood]->Fill(dy);
413 fpDPtS[kGood]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
414 }
415
416 } // good matches
417 } // found possible match
418 } // upper
419 } // tracks 1
420 PostData(1, fHists);
421}
422
423//________________________________________________________________________
424void AliAnalysisTaskCosmic::Terminate(Option_t *)
425{
426}