]>
Commit | Line | Data |
---|---|---|
16a89872 | 1 | // ##### TO RUN THIS MACRO: |
2 | // bash$ aliroot (due to AliInfo, AliError, ...) | |
c52c2132 | 3 | // root[0] gSystem->Load("libANALYSIS"); |
16a89872 | 4 | // IN case you do not have the include path set to AliRoot includes: |
5 | // root[1] gSystem->AddIncludePath("-I\"$ALICE_ROOT/include\""); | |
6 | // root[2] .L testEvent.C+; | |
7 | // root[3] generate(); | |
8 | // root[4] filter_reco(); | |
9 | ||
10 | #include "TClonesArray.h" | |
11 | #include "TChain.h" | |
12 | #include "TH1.h" | |
c52c2132 | 13 | #include "TMath.h" |
16a89872 | 14 | #include "TCanvas.h" |
15 | #include "testEvent.h" | |
16 | ||
17 | //============= First step: generate events | |
18 | void generate() | |
19 | { | |
20 | // Simple event generation | |
c52c2132 | 21 | AliAnalysisManager *mgr = new AliAnalysisManager("generate"); |
16a89872 | 22 | TaskGenerate *task = new TaskGenerate("gener"); |
23 | mgr->AddTask(task); | |
24 | ||
25 | if (mgr->InitAnalysis()) { | |
26 | mgr->PrintStatus(); | |
c52c2132 | 27 | mgr->StartAnalysis(); |
16a89872 | 28 | } |
29 | delete mgr; | |
30 | } | |
31 | ||
32 | //============= Second step: filter gammas; use TSelector functionality | |
33 | void filter_reco() | |
34 | { | |
35 | // Filter the input events having more than 100 gammas. Reconstruct pi0's | |
36 | // From gammas coming from the same vertex. | |
37 | // Get the input data as a chain | |
38 | TChain *chain = new TChain("T"); | |
39 | chain->Add("event02000.root"); | |
40 | chain->Add("event04000.root"); | |
41 | chain->Add("event06000.root"); | |
42 | chain->Add("event08000.root"); | |
43 | chain->Add("event10000.root"); | |
44 | // Create an analysis manager | |
c52c2132 | 45 | AliAnalysisManager *mgr = new AliAnalysisManager("testEvent"); |
16a89872 | 46 | // Create a filter task and register it |
47 | TaskFilter *task1 = new TaskFilter("TaskFilter"); | |
48 | mgr->AddTask(task1); | |
49 | // Create a reco task and register it | |
c52c2132 | 50 | |
16a89872 | 51 | TaskRecoPi0 *task2 = new TaskRecoPi0("TaskRecoPi0"); |
52 | mgr->AddTask(task2); | |
53 | // Create containers for input/output | |
54 | AliAnalysisDataContainer *cinput = mgr->CreateContainer("input0", | |
55 | TTree::Class(), AliAnalysisManager::kInputContainer); | |
56 | AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("output0", | |
57 | TTree::Class(), AliAnalysisManager::kOutputContainer); | |
58 | AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("output1", | |
59 | TList::Class(), AliAnalysisManager::kOutputContainer); | |
60 | AliAnalysisDataContainer *coutput = mgr->CreateContainer("output", | |
c52c2132 | 61 | TH1::Class(), AliAnalysisManager::kOutputContainer, "output.root"); |
16a89872 | 62 | // Connect containers to the task input/outputs |
63 | mgr->ConnectInput(task1,0,cinput); | |
64 | mgr->ConnectOutput(task1,0,coutput1); | |
65 | mgr->ConnectOutput(task1,1,coutput2); | |
66 | mgr->ConnectInput(task2,0,coutput1); | |
67 | mgr->ConnectOutput(task2,0,coutput); | |
c52c2132 | 68 | |
16a89872 | 69 | // Connect input data |
c52c2132 | 70 | // cinput->SetData(chain); |
16a89872 | 71 | // Init analysis and start event loop |
72 | if (mgr->InitAnalysis()) { | |
73 | mgr->PrintStatus(); | |
c52c2132 | 74 | mgr->StartAnalysis("local",chain); |
16a89872 | 75 | } |
c52c2132 | 76 | delete mgr; |
16a89872 | 77 | } |
78 | ||
79 | ClassImp(TaskGenerate) | |
80 | ||
81 | //______________________________________________________________________________ | |
82 | void TaskGenerate::Exec(Option_t *) | |
83 | { | |
84 | // Generate 10k events in 5 files | |
85 | AnaEvent::CreateEvents(2000, "event02000.root"); | |
86 | AnaEvent::CreateEvents(2000, "event04000.root"); | |
87 | AnaEvent::CreateEvents(2000, "event06000.root"); | |
88 | AnaEvent::CreateEvents(2000, "event08000.root"); | |
89 | AnaEvent::CreateEvents(2000, "event10000.root"); | |
90 | } | |
91 | ||
92 | ClassImp(TaskFilter) | |
93 | ||
94 | //______________________________________________________________________________ | |
95 | TaskFilter::TaskFilter(const char *name) | |
96 | :AliAnalysisTask(name,""), fEvent(0), fOutput(0), fList(0), fHist1(0), fHist2(0) | |
97 | { | |
98 | // Constructor | |
99 | // Input slot #0 works with a tree of events | |
100 | DefineInput(0, TTree::Class()); | |
101 | ||
102 | // Output slot #0 writes into a tree, #1 produces a hisogram | |
103 | DefineOutput(0, TTree::Class()); | |
104 | DefineOutput(1, TList::Class()); | |
105 | } | |
106 | ||
107 | //______________________________________________________________________________ | |
c52c2132 | 108 | void TaskFilter::ConnectInputData(Option_t *) |
16a89872 | 109 | { |
110 | // Initialize branches. | |
c52c2132 | 111 | printf(" ConnectInputData of task %s\n", GetName()); |
112 | char ** address = (char **)GetBranchAddress(0, "event"); | |
113 | if (address) { | |
16a89872 | 114 | // One should first check if the branch address was taken by some other task |
c52c2132 | 115 | fEvent = (AnaEvent*)(*address); |
116 | } else { | |
117 | fEvent = new AnaEvent(); | |
118 | SetBranchAddress(0, "event", &fEvent); | |
119 | } | |
120 | } | |
121 | ||
122 | //______________________________________________________________________________ | |
123 | void TaskFilter::CreateOutputObjects() | |
124 | { | |
125 | printf(" CreateOutputObjects of task %s\n", GetName()); | |
126 | if (!fList) { | |
16a89872 | 127 | fOutput = new TTree("TGAM", "gammas"); |
128 | TBranch *branch = fOutput->Branch("event", &fEvent, 18000,1); | |
129 | branch->SetAutoDelete(kFALSE); | |
16a89872 | 130 | fList = new TList(); |
131 | fHist1 = new TH1I("ntracks", "Number of tracks per event", 100, 0, 1000); | |
132 | fHist1->SetLineColor(kRed); | |
133 | fHist2 = new TH1I("ngammas", "Number of gammas per event", 100, 0, 1000); | |
134 | fHist2->SetLineColor(kBlue); | |
135 | fList->Add(fHist1); | |
136 | fList->Add(fHist2); | |
137 | } | |
138 | } | |
139 | ||
140 | //______________________________________________________________________________ | |
141 | void TaskFilter::Exec(Option_t *) | |
142 | { | |
143 | // Filtering. | |
144 | TTree *tinput = (TTree*)GetInputData(0); | |
145 | Long64_t ientry = tinput->GetReadEntry(); | |
146 | // In this case fEvent address is already connected to the input | |
147 | if (!fEvent) return; | |
148 | // First check multiplicity | |
149 | Int_t ntracks = fEvent->GetNtracks(); | |
150 | // Loop tracks and get rid of non-gammas | |
151 | AnaTrack *track; | |
152 | Int_t igamma = 0; | |
153 | for (Int_t i=0; i<ntracks; i++) { | |
154 | track = fEvent->GetTrack(i); | |
155 | if (track->GetMass() < 1.e-3) igamma++; | |
156 | } | |
157 | if (ientry%100 == 0) printf("TaskFilter -> Event %lld: %i tracks filtered %i gammas\n", ientry, ntracks, igamma); | |
158 | fHist1->Fill(ntracks); | |
159 | fHist2->Fill(igamma); | |
160 | if (igamma > 100) { | |
161 | fOutput->Fill(); | |
162 | PostData(0, fOutput); | |
163 | } | |
164 | PostData(1, fList); | |
165 | } | |
166 | ||
167 | //______________________________________________________________________________ | |
168 | void TaskFilter::Terminate(Option_t *) | |
169 | { | |
170 | // Draw some histogram at the end. | |
171 | if (!gROOT->IsBatch()) { | |
172 | fHist1->SetMaximum(2500); | |
173 | fHist1->Draw(); | |
174 | fHist2->Draw("SAME"); | |
175 | } | |
176 | } | |
177 | ||
178 | ClassImp(TaskRecoPi0) | |
179 | ||
180 | //______________________________________________________________________________ | |
181 | TaskRecoPi0::TaskRecoPi0(const char *name) | |
182 | :AliAnalysisTask(name,""), fEvent(0), fGammas(0), fPions(0), fHist(0) | |
183 | { | |
184 | // Constructor | |
185 | // Input slot #0 works with a tree of events | |
186 | DefineInput(0, TTree::Class()); | |
187 | ||
188 | // Output slot #1 produces a hisogram | |
189 | DefineOutput(0, TH1::Class()); | |
190 | } | |
191 | ||
192 | //______________________________________________________________________________ | |
193 | TaskRecoPi0::~TaskRecoPi0() | |
194 | { | |
195 | // Dtor. | |
196 | if (fEvent) delete fEvent; | |
197 | if (fGammas) delete fGammas; | |
198 | if (fPions) delete fPions; | |
199 | } | |
200 | ||
201 | //______________________________________________________________________________ | |
c52c2132 | 202 | void TaskRecoPi0::ConnectInputData(Option_t *) |
16a89872 | 203 | { |
204 | // Initialize branches. | |
c52c2132 | 205 | printf(" ConnectInputData for task %s\n", GetName()); |
16a89872 | 206 | if (!fEvent) { |
207 | // One should first check if the branch address was taken by some other task | |
208 | char ** address = (char **)GetBranchAddress(0, "event"); | |
209 | if (address) fEvent = (AnaEvent*)(*address); | |
210 | if (!fEvent) { | |
211 | fEvent = new AnaEvent(); | |
212 | SetBranchAddress(0, "event", &fEvent); | |
213 | } | |
16a89872 | 214 | } |
c52c2132 | 215 | } |
216 | ||
217 | //______________________________________________________________________________ | |
218 | void TaskRecoPi0::CreateOutputObjects() | |
219 | { | |
220 | printf(" CreateOutputObjects of task %s\n", GetName()); | |
16a89872 | 221 | if (!fHist) { |
c52c2132 | 222 | fGammas = new TObjArray(); |
223 | fPions = new TObjArray(); | |
16a89872 | 224 | fHist = new TH1F("Pt_pi0", "Pt distribution for pi0's", 100, 0., 10.); |
225 | fHist->SetLineColor(kRed); | |
226 | } | |
227 | } | |
228 | ||
229 | //______________________________________________________________________________ | |
230 | void TaskRecoPi0::Exec(Option_t *) | |
231 | { | |
232 | // Reconstruct Pi0's for one event | |
233 | AnaTrack *track = 0; | |
234 | Int_t ntracks = fEvent->GetNtracks(); | |
235 | Int_t ngamma = 0; | |
236 | Int_t i,j; | |
237 | // Clear containers | |
238 | fGammas->Clear(); | |
239 | fPions->Delete(); | |
240 | fPions->Clear(); | |
241 | // Loop tracks and move gammas to gamma container | |
242 | for (i=0; i<ntracks; i++) { | |
243 | track = fEvent->GetTrack(i); | |
244 | if (track->GetMass() < 1.e-3) { | |
245 | fGammas->Add(track); | |
246 | ngamma++; | |
247 | } | |
248 | } | |
249 | printf("TaskRecoPi0 -> Tracks %i \n", ntracks); | |
250 | ||
251 | // Loop gammas and check vertex position | |
252 | Double_t v1[3], v2[3]; | |
253 | AnaTrack *tracko = 0; | |
254 | Double_t cutoff = 0.001; | |
255 | for (i=0; i<ngamma-1; i++) { | |
256 | track = (AnaTrack*)fGammas->At(i); | |
257 | v1[0] = track->GetVertex(0); | |
258 | v1[1] = track->GetVertex(1); | |
259 | v1[2] = track->GetVertex(2); | |
260 | for (j=i+1; j<ngamma; j++) { | |
261 | tracko = (AnaTrack*)fGammas->At(j); | |
262 | v2[0] = tracko->GetVertex(0); | |
263 | v2[1] = tracko->GetVertex(1); | |
264 | v2[2] = tracko->GetVertex(2); | |
265 | Double_t dist2 = (v2[0]-v1[0])*(v2[0]-v1[0])+ | |
266 | (v2[1]-v1[1])*(v2[1]-v1[1])+ | |
267 | (v2[2]-v1[2])*(v2[2]-v1[2]); | |
268 | if (dist2>cutoff*cutoff) continue; | |
269 | // We have found a pair candidate | |
270 | Double_t px = track->GetPx()+tracko->GetPx(); | |
271 | Double_t py = track->GetPy()+tracko->GetPy(); | |
272 | Double_t pz = track->GetPz()+tracko->GetPz(); | |
273 | track = new AnaTrack(px,py,pz,0.135,0,v1[0],v1[1],v1[2]); | |
274 | fPions->Add(track); | |
275 | fHist->Fill(track->GetPt()); | |
276 | break; | |
277 | } | |
278 | } | |
279 | PostData(0,fHist); | |
280 | } | |
281 | ||
282 | //______________________________________________________________________________ | |
283 | void TaskRecoPi0::Terminate(Option_t *) | |
284 | { | |
285 | // Draw some histogram at the end. | |
286 | if (!gROOT->IsBatch()) { | |
287 | new TCanvas("pi0", "Pt for pi0's", 800,600); | |
288 | fHist->Draw(); | |
289 | } | |
290 | } | |
291 | ||
292 | ||
293 | ClassImp(AnaTrack) | |
294 | ||
295 | //______________________________________________________________________________ | |
296 | AnaTrack::AnaTrack(Double_t random, Double_t *vertex) | |
297 | { | |
298 | // Constructor | |
299 | Int_t itype; | |
300 | if (random<0.3) itype = 1; // pi+ | |
301 | else if (random<0.6) itype = 2; // pi- | |
302 | else if (random<0.9) itype = 3; // p | |
303 | else if (random<0.95) itype = 4; // pi0 | |
304 | else itype = 5; // gamma | |
305 | gRandom->Rannor(fPx, fPy); | |
306 | Double_t vert_width = 0.1; | |
307 | ||
308 | switch (itype) { | |
309 | case 1: | |
310 | fMass = 0.13957 + gRandom->Gaus(0.,0.001); | |
311 | fCharge = 1; | |
312 | break; | |
313 | case 2: | |
314 | fMass = 0.13957 + gRandom->Gaus(0.,0.001); | |
315 | fCharge = -1; | |
316 | break; | |
317 | case 3: | |
318 | fMass = 0.938 + gRandom->Gaus(0.,0.002); | |
319 | fCharge = 1; | |
320 | fPx *= 0.15; | |
321 | fPy *= 0.15; | |
322 | break; | |
323 | case 4: | |
324 | fMass = 0.135 + gRandom->Gaus(0.,0.001); | |
325 | fCharge = 0; | |
326 | fPx *= 0.8; | |
327 | fPy *= 0.8; | |
328 | vert_width = 10.; | |
329 | break; | |
330 | case 5: | |
331 | fMass = 0.; | |
332 | fCharge = 0; | |
333 | fPx *= 0.5; | |
334 | fPy *= 1.5; | |
335 | break; | |
336 | }; | |
337 | fPz = gRandom->Gaus(4., 2.); | |
338 | if (vertex) { | |
339 | fVertex[0] = vertex[0]; | |
340 | fVertex[1] = vertex[1]; | |
341 | fVertex[2] = vertex[2]; | |
342 | } else { | |
343 | fVertex[0] = gRandom->Gaus(0,vert_width); | |
344 | fVertex[1] = gRandom->Gaus(0,vert_width); | |
345 | fVertex[2] = gRandom->Gaus(0,0.01); | |
346 | } | |
347 | } | |
348 | ||
349 | //______________________________________________________________________________ | |
350 | Bool_t AnaTrack::Decay(Double_t &px1, Double_t &py1, Double_t &pz1, | |
351 | Double_t &px2, Double_t &py2, Double_t &pz2) | |
352 | { | |
353 | // Decay a pi0 in 2 gammas. | |
354 | if (fCharge != 0) return kFALSE; | |
355 | if (fMass<0.132 || fMass>0.138) return kFALSE; | |
356 | Double_t phi1 = 2.*TMath::Pi()*gRandom->Rndm(); // [0,2*pi] | |
357 | Double_t phi2 = phi1 + TMath::Pi(); | |
358 | if (phi2 > 2.*TMath::Pi()) phi2 -= 2.*TMath::Pi(); | |
359 | Double_t r2 = gRandom->Rndm(); | |
360 | Double_t theta1 = TMath::ACos(1.-2.*r2); | |
361 | Double_t p0 = GetP(); | |
362 | Double_t m0 = 0.135; | |
363 | Double_t p1 = 0.5*(p0*(1-2*r2)+TMath::Sqrt(p0*p0*(1-2*r2)*(1-2*r2)+2*m0*m0)); | |
364 | Double_t p2 = TMath::Sqrt(p0*p0+m0*m0-p1*p1); | |
365 | Double_t theta2 = TMath::ACos((p0-p1*(1-2*r2))/p2); | |
366 | // Px, Py and Pz in the reference frame of the pion | |
367 | px1 = p1 * TMath::Sin(theta1)*TMath::Cos(phi1); | |
368 | px2 = p2 * TMath::Sin(theta2)*TMath::Cos(phi2); | |
369 | py1 = p1 * TMath::Sin(theta1)*TMath::Sin(phi1); | |
370 | py2 = p2 * TMath::Sin(theta2)*TMath::Sin(phi2); | |
371 | pz1 = p1 * TMath::Cos(theta1); | |
372 | pz2 = p2 * TMath::Cos(theta2); | |
373 | Double_t phi = TMath::ATan2(fPy, fPx) * TMath::RadToDeg(); | |
374 | Double_t theta = TMath::ACos(GetPt()/p0) * TMath::RadToDeg(); | |
375 | TGeoRotation r("rot", phi+90., theta, 0.); | |
376 | Double_t loc[3], vect[3]; | |
377 | p0 = TMath::Sqrt(px1*px1+py1*py1+pz1*pz1); | |
378 | loc[0] = px1/p0; | |
379 | loc[1] = py1/p0; | |
380 | loc[2] = pz1/p0; | |
381 | r.LocalToMasterVect(loc, vect); | |
382 | px1 = vect[0]*p0; | |
383 | py1 = vect[1]*p0; | |
384 | pz1 = vect[2]*p0; | |
385 | // t1 = new AnaTrack(1., fVertex); | |
386 | p0 = TMath::Sqrt(px2*px2+py2*py2+pz2*pz2); | |
387 | loc[0] = px2/p0; | |
388 | loc[1] = py2/p0; | |
389 | loc[2] = pz2/p0; | |
390 | r.LocalToMasterVect(loc, vect); | |
391 | px2 = vect[0]*p0; | |
392 | py2 = vect[1]*p0; | |
393 | pz2 = vect[2]*p0; | |
394 | // t2 = new AnaTrack(1., fVertex); | |
395 | return kTRUE; | |
396 | } | |
397 | ||
398 | //______________________________________________________________________________ | |
399 | ClassImp(AnaEvent) | |
400 | ||
401 | TClonesArray *AnaEvent::fgTracks = 0; | |
402 | ||
403 | //______________________________________________________________________________ | |
404 | AnaEvent::AnaEvent() | |
405 | { | |
406 | // Ctor | |
407 | if (!fgTracks) fgTracks = new TClonesArray("AnaTrack", 2000); | |
408 | fTracks = fgTracks; | |
409 | fEventNumber = 0; | |
410 | fNtracks = 0; | |
411 | } | |
412 | ||
413 | //______________________________________________________________________________ | |
414 | AnaTrack *AnaEvent::AddTrack(Double_t rnd, Double_t *vert) | |
415 | { | |
416 | // Add a random track | |
417 | // printf("track %d\n", fNtracks); | |
418 | TClonesArray &tracks = *fTracks; | |
419 | AnaTrack *track = new(tracks[fNtracks++]) AnaTrack(rnd, vert); | |
420 | return track; | |
421 | } | |
422 | ||
423 | //______________________________________________________________________________ | |
424 | void AnaEvent::Clear(Option_t *) | |
425 | { | |
426 | // Clears current event. | |
427 | fTracks->Clear(); | |
428 | fNtracks = 0; | |
429 | fEventNumber = 0; | |
430 | } | |
431 | ||
432 | //______________________________________________________________________________ | |
433 | Int_t AnaEvent::Build(Int_t ev) | |
434 | { | |
435 | // Create a random event | |
436 | Clear(); | |
437 | fEventNumber = ev; | |
438 | Int_t ntracks = Int_t(gRandom->Gaus(500., 100.)); | |
439 | if (ntracks < 1) ntracks = 1; | |
440 | if (ntracks>1000) ntracks = 1000; | |
441 | AnaTrack *track, *track0; | |
442 | for (Int_t i=0; i<ntracks; i++) { | |
443 | Double_t rnd = gRandom->Rndm(); | |
444 | if (rnd>=0.90 && rnd<0.95) { | |
445 | // Pi0 -> decay the track in 2 gammas | |
446 | track0 = new AnaTrack(rnd); | |
447 | Double_t vert[3]; | |
448 | vert[0] = track0->GetVertex(0); | |
449 | vert[1] = track0->GetVertex(1); | |
450 | vert[2] = track0->GetVertex(2); | |
451 | Double_t px1,py1,pz1,px2,py2,pz2; | |
452 | if (track0->Decay(px1,py1,pz1,px2,py2,pz2)) { | |
453 | track = AddTrack(1.,vert); | |
454 | track->SetPx(px1); | |
455 | track->SetPy(py1); | |
456 | track->SetPz(pz1); | |
457 | track = AddTrack(1.,vert); | |
458 | track->SetPx(px2); | |
459 | track->SetPy(py2); | |
460 | track->SetPz(pz2); | |
461 | } | |
462 | delete track0; | |
463 | } else { | |
464 | track = AddTrack(rnd); | |
465 | } | |
466 | } | |
467 | return fNtracks; | |
468 | } | |
469 | ||
470 | //______________________________________________________________________________ | |
471 | void AnaEvent::CreateEvents(Int_t nevents, const char *filename) | |
472 | { | |
473 | // Create nevents in one tree and put them in filename. | |
474 | TFile *hfile = new TFile(filename, "RECREATE", "Some AnaEvents..."); | |
475 | TTree *tree = new TTree("T", "Tree of AnaEvents"); | |
476 | tree->SetAutoSave(1000000000); // autosave when 1 Gbyte written | |
477 | Int_t bufsize = 16000; | |
478 | AnaEvent *event = new AnaEvent(); | |
479 | TBranch *branch = tree->Branch("event", &event, bufsize,1); | |
480 | branch->SetAutoDelete(kFALSE); | |
481 | ||
482 | for (Int_t ev=0; ev<nevents; ev++) { | |
483 | Int_t ntracks = event->Build(ev); | |
484 | if (ev%100 == 0) printf("event: %d ntracks=%d\n", ev, ntracks); | |
485 | tree->Fill(); | |
486 | } | |
487 | hfile->Write(); | |
488 | tree->Print(); | |
489 | hfile->Close(); | |
490 | } | |
491 |