coverity
[u/mrichter/AliRoot.git] / ANALYSIS / testEvent.C
CommitLineData
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
18void 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
33void 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
79ClassImp(TaskGenerate)
80
81//______________________________________________________________________________
82void 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
92ClassImp(TaskFilter)
93
94//______________________________________________________________________________
95TaskFilter::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 108void 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//______________________________________________________________________________
123void 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//______________________________________________________________________________
141void 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//______________________________________________________________________________
168void 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
178ClassImp(TaskRecoPi0)
179
180//______________________________________________________________________________
181TaskRecoPi0::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//______________________________________________________________________________
193TaskRecoPi0::~TaskRecoPi0()
194{
195// Dtor.
196 if (fEvent) delete fEvent;
197 if (fGammas) delete fGammas;
198 if (fPions) delete fPions;
199}
200
201//______________________________________________________________________________
c52c2132 202void 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//______________________________________________________________________________
218void 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//______________________________________________________________________________
230void 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//______________________________________________________________________________
283void 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
293ClassImp(AnaTrack)
294
295//______________________________________________________________________________
296AnaTrack::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//______________________________________________________________________________
350Bool_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//______________________________________________________________________________
399ClassImp(AnaEvent)
400
401TClonesArray *AnaEvent::fgTracks = 0;
402
403//______________________________________________________________________________
404AnaEvent::AnaEvent()
405{
406// Ctor
407 if (!fgTracks) fgTracks = new TClonesArray("AnaTrack", 2000);
408 fTracks = fgTracks;
409 fEventNumber = 0;
410 fNtracks = 0;
411}
412
413//______________________________________________________________________________
414AnaTrack *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//______________________________________________________________________________
424void AnaEvent::Clear(Option_t *)
425{
426// Clears current event.
427 fTracks->Clear();
428 fNtracks = 0;
429 fEventNumber = 0;
430}
431
432//______________________________________________________________________________
433Int_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//______________________________________________________________________________
471void 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