1 // ##### TO RUN THIS MACRO:
2 // bash$ aliroot (due to AliInfo, AliError, ...)
3 // root[0] gSystem->Load("libANALYSIS");
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+;
8 // root[4] filter_reco();
10 #include "TClonesArray.h"
15 #include "testEvent.h"
17 //============= First step: generate events
20 // Simple event generation
21 AliAnalysisManager *mgr = new AliAnalysisManager("generate");
22 TaskGenerate *task = new TaskGenerate("gener");
25 if (mgr->InitAnalysis()) {
32 //============= Second step: filter gammas; use TSelector functionality
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
45 AliAnalysisManager *mgr = new AliAnalysisManager("testEvent");
46 // Create a filter task and register it
47 TaskFilter *task1 = new TaskFilter("TaskFilter");
49 // Create a reco task and register it
51 TaskRecoPi0 *task2 = new TaskRecoPi0("TaskRecoPi0");
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",
61 TH1::Class(), AliAnalysisManager::kOutputContainer, "output.root");
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);
70 // cinput->SetData(chain);
71 // Init analysis and start event loop
72 if (mgr->InitAnalysis()) {
74 mgr->StartAnalysis("local",chain);
79 ClassImp(TaskGenerate)
81 //______________________________________________________________________________
82 void TaskGenerate::Exec(Option_t *)
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");
94 //______________________________________________________________________________
95 TaskFilter::TaskFilter(const char *name)
96 :AliAnalysisTask(name,""), fEvent(0), fOutput(0), fList(0), fHist1(0), fHist2(0)
99 // Input slot #0 works with a tree of events
100 DefineInput(0, TTree::Class());
102 // Output slot #0 writes into a tree, #1 produces a hisogram
103 DefineOutput(0, TTree::Class());
104 DefineOutput(1, TList::Class());
107 //______________________________________________________________________________
108 void TaskFilter::ConnectInputData(Option_t *)
110 // Initialize branches.
111 printf(" ConnectInputData of task %s\n", GetName());
112 char ** address = (char **)GetBranchAddress(0, "event");
114 // One should first check if the branch address was taken by some other task
115 fEvent = (AnaEvent*)(*address);
117 fEvent = new AnaEvent();
118 SetBranchAddress(0, "event", &fEvent);
122 //______________________________________________________________________________
123 void TaskFilter::CreateOutputObjects()
125 printf(" CreateOutputObjects of task %s\n", GetName());
127 fOutput = new TTree("TGAM", "gammas");
128 TBranch *branch = fOutput->Branch("event", &fEvent, 18000,1);
129 branch->SetAutoDelete(kFALSE);
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);
140 //______________________________________________________________________________
141 void TaskFilter::Exec(Option_t *)
144 TTree *tinput = (TTree*)GetInputData(0);
145 Long64_t ientry = tinput->GetReadEntry();
146 // In this case fEvent address is already connected to the input
148 // First check multiplicity
149 Int_t ntracks = fEvent->GetNtracks();
150 // Loop tracks and get rid of non-gammas
153 for (Int_t i=0; i<ntracks; i++) {
154 track = fEvent->GetTrack(i);
155 if (track->GetMass() < 1.e-3) igamma++;
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);
162 PostData(0, fOutput);
167 //______________________________________________________________________________
168 void TaskFilter::Terminate(Option_t *)
170 // Draw some histogram at the end.
171 if (!gROOT->IsBatch()) {
172 fHist1->SetMaximum(2500);
174 fHist2->Draw("SAME");
178 ClassImp(TaskRecoPi0)
180 //______________________________________________________________________________
181 TaskRecoPi0::TaskRecoPi0(const char *name)
182 :AliAnalysisTask(name,""), fEvent(0), fGammas(0), fPions(0), fHist(0)
185 // Input slot #0 works with a tree of events
186 DefineInput(0, TTree::Class());
188 // Output slot #1 produces a hisogram
189 DefineOutput(0, TH1::Class());
192 //______________________________________________________________________________
193 TaskRecoPi0::~TaskRecoPi0()
196 if (fEvent) delete fEvent;
197 if (fGammas) delete fGammas;
198 if (fPions) delete fPions;
201 //______________________________________________________________________________
202 void TaskRecoPi0::ConnectInputData(Option_t *)
204 // Initialize branches.
205 printf(" ConnectInputData for task %s\n", GetName());
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);
211 fEvent = new AnaEvent();
212 SetBranchAddress(0, "event", &fEvent);
217 //______________________________________________________________________________
218 void TaskRecoPi0::CreateOutputObjects()
220 printf(" CreateOutputObjects of task %s\n", GetName());
222 fGammas = new TObjArray();
223 fPions = new TObjArray();
224 fHist = new TH1F("Pt_pi0", "Pt distribution for pi0's", 100, 0., 10.);
225 fHist->SetLineColor(kRed);
229 //______________________________________________________________________________
230 void TaskRecoPi0::Exec(Option_t *)
232 // Reconstruct Pi0's for one event
234 Int_t ntracks = fEvent->GetNtracks();
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) {
249 printf("TaskRecoPi0 -> Tracks %i \n", ntracks);
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]);
275 fHist->Fill(track->GetPt());
282 //______________________________________________________________________________
283 void TaskRecoPi0::Terminate(Option_t *)
285 // Draw some histogram at the end.
286 if (!gROOT->IsBatch()) {
287 new TCanvas("pi0", "Pt for pi0's", 800,600);
295 //______________________________________________________________________________
296 AnaTrack::AnaTrack(Double_t random, Double_t *vertex)
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;
310 fMass = 0.13957 + gRandom->Gaus(0.,0.001);
314 fMass = 0.13957 + gRandom->Gaus(0.,0.001);
318 fMass = 0.938 + gRandom->Gaus(0.,0.002);
324 fMass = 0.135 + gRandom->Gaus(0.,0.001);
337 fPz = gRandom->Gaus(4., 2.);
339 fVertex[0] = vertex[0];
340 fVertex[1] = vertex[1];
341 fVertex[2] = vertex[2];
343 fVertex[0] = gRandom->Gaus(0,vert_width);
344 fVertex[1] = gRandom->Gaus(0,vert_width);
345 fVertex[2] = gRandom->Gaus(0,0.01);
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)
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();
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);
381 r.LocalToMasterVect(loc, vect);
385 // t1 = new AnaTrack(1., fVertex);
386 p0 = TMath::Sqrt(px2*px2+py2*py2+pz2*pz2);
390 r.LocalToMasterVect(loc, vect);
394 // t2 = new AnaTrack(1., fVertex);
398 //______________________________________________________________________________
401 TClonesArray *AnaEvent::fgTracks = 0;
403 //______________________________________________________________________________
407 if (!fgTracks) fgTracks = new TClonesArray("AnaTrack", 2000);
413 //______________________________________________________________________________
414 AnaTrack *AnaEvent::AddTrack(Double_t rnd, Double_t *vert)
416 // Add a random track
417 // printf("track %d\n", fNtracks);
418 TClonesArray &tracks = *fTracks;
419 AnaTrack *track = new(tracks[fNtracks++]) AnaTrack(rnd, vert);
423 //______________________________________________________________________________
424 void AnaEvent::Clear(Option_t *)
426 // Clears current event.
432 //______________________________________________________________________________
433 Int_t AnaEvent::Build(Int_t ev)
435 // Create a random event
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);
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);
457 track = AddTrack(1.,vert);
464 track = AddTrack(rnd);
470 //______________________________________________________________________________
471 void AnaEvent::CreateEvents(Int_t nevents, const char *filename)
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);
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);