fixed compiler problem
[u/mrichter/AliRoot.git] / ANALYSIS / testEvent.C
CommitLineData
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
17void 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
32void 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
76ClassImp(TaskGenerate)
77
78//______________________________________________________________________________
79void 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
89ClassImp(TaskFilter)
90
91//______________________________________________________________________________
92TaskFilter::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//______________________________________________________________________________
105void 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//______________________________________________________________________________
138void 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//______________________________________________________________________________
165void 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
175ClassImp(TaskRecoPi0)
176
177//______________________________________________________________________________
178TaskRecoPi0::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//______________________________________________________________________________
190TaskRecoPi0::~TaskRecoPi0()
191{
192// Dtor.
193 if (fEvent) delete fEvent;
194 if (fGammas) delete fGammas;
195 if (fPions) delete fPions;
196}
197
198//______________________________________________________________________________
199void 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//______________________________________________________________________________
221void 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//______________________________________________________________________________
274void 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
284ClassImp(AnaTrack)
285
286//______________________________________________________________________________
287AnaTrack::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//______________________________________________________________________________
341Bool_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//______________________________________________________________________________
390ClassImp(AnaEvent)
391
392TClonesArray *AnaEvent::fgTracks = 0;
393
394//______________________________________________________________________________
395AnaEvent::AnaEvent()
396{
397// Ctor
398 if (!fgTracks) fgTracks = new TClonesArray("AnaTrack", 2000);
399 fTracks = fgTracks;
400 fEventNumber = 0;
401 fNtracks = 0;
402}
403
404//______________________________________________________________________________
405AnaTrack *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//______________________________________________________________________________
415void AnaEvent::Clear(Option_t *)
416{
417// Clears current event.
418 fTracks->Clear();
419 fNtracks = 0;
420 fEventNumber = 0;
421}
422
423//______________________________________________________________________________
424Int_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//______________________________________________________________________________
462void 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