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