1 Int_t diaglevel=2; // 1->Hits, 2->Spectra, 3->Statistics
4 void RICHpadtestC (Int_t evNumber1=0,Int_t evNumber2=0)
6 /////////////////////////////////////////////////////////////////////////
7 // This macro is a small example of a ROOT macro
8 // illustrating how to read the output of GALICE
9 // and do some analysis.
11 /////////////////////////////////////////////////////////////////////////
14 Int_t NpadX = 162; // number of pads on X
15 Int_t NpadY = 162; // number of pads on Y
18 for (Int_t i=0;i<NpadX;i++) {
19 for (Int_t j=0;j<NpadY;j++) {
26 // Dynamically link some shared libs
28 if (gClassTable->GetID("AliRun") < 0) {
29 gROOT->LoadMacro("loadlibs.C");
39 // Connect the Root Galice file containing Geometry, Kine and Hits
41 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
42 if (!file) file = new TFile("galice.root","UPDATE");
44 // Get AliRun object from file or create it if not on file
47 gAlice = (AliRun*)file->Get("gAlice");
48 if (gAlice) printf("AliRun object found on file\n");
49 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
53 gAlice = (AliRun*)file->Get("gAlice");
54 if (gAlice) printf("AliRun object found on file\n");
55 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
58 // Create some histograms
65 TH1F *pionspectra1 = new TH1F("pionspectra1","Pion Spectra",200,-4,2);
66 TH1F *pionspectra2 = new TH1F("pionspectra2","Pion Spectra",200,-4,2);
67 TH1F *pionspectra3 = new TH1F("pionspectra3","Pion Spectra",200,-4,2);
68 TH1F *protonspectra1 = new TH1F("protonspectra1","Proton Spectra",200,-4,2);
69 TH1F *protonspectra2 = new TH1F("protonspectra2","Proton Spectra",200,-4,2);
70 TH1F *protonspectra3 = new TH1F("protonspectra3","Proton Spectra",200,-4,2);
71 TH1F *kaonspectra1 = new TH1F("kaonspectra1","Kaon Spectra",100,-4,2);
72 TH1F *kaonspectra2 = new TH1F("kaonspectra2","Kaon Spectra",100,-4,2);
73 TH1F *kaonspectra3 = new TH1F("kaonspectra3","Kaon Spectra",100,-4,2);
74 TH1F *electronspectra1 = new TH1F("electronspectra1","Electron Spectra",100,-4,2);
75 TH1F *electronspectra2 = new TH1F("electronspectra2","Electron Spectra",100,-4,2);
76 TH1F *electronspectra3 = new TH1F("electronspectra3","Electron Spectra",100,-4,2);
77 TH1F *muonspectra1 = new TH1F("muonspectra1","Muon Spectra",100,-4,2);
78 TH1F *muonspectra2 = new TH1F("muonspectra2","Muon Spectra",100,-4,2);
79 TH1F *muonspectra3 = new TH1F("muonspectra3","Muon Spectra",100,-4,2);
80 TH1F *neutronspectra1 = new TH1F("neutronspectra1","Neutron Spectra",100,-4,2);
81 TH1F *neutronspectra2 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
82 TH1F *neutronspectra3 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
83 TH1F *chargedspectra1 = new TH1F("chargedspectra1","Charged particles above 1 GeV Spectra",100,-1,3);
84 TH1F *chargedspectra2 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
85 TH1F *chargedspectra3 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
86 TH2F *production = new TH2F("production","Mother production vertices",100,-300,300,100,0,600);
91 // Start loop over events
101 Int_t pion=0, kaon=0, proton=0, electron=0, positron=0, neutron=0, highneutrons=0, muon=0;
102 Int_t chargedpions=0,primarypions=0,highprimarypions=0,chargedkaons=0,primarykaons=0,highprimarykaons=0;
103 Int_t chargedmuons=0, photons=0, primaryphotons=0, highprimaryphotons=0;
107 for (int nev=0; nev<= evNumber2; nev++) {
108 Int_t nparticles = gAlice->GetEvent(nev);
111 printf ("Event number : %d\n",nev);
112 printf ("Number of particles: %d\n",nparticles);
113 if (nev < evNumber1) continue;
114 if (nparticles <= 0) return;
116 // Get pointers to RICH detector and Hits containers
118 AliRICH *RICH = (AliRICH*)gAlice->GetDetector("RICH");
119 Int_t nent=(Int_t)gAlice->TreeR()->GetEntries();
120 gAlice->TreeR()->GetEvent(nent-1);
121 TClonesArray *Rawclusters = RICH->RawClustAddress(2); // Raw clusters branch
122 Int_t nrawclusters = Rawclusters->GetEntriesFast();
123 TTree *TH = gAlice->TreeH();
124 Int_t ntracks = TH->GetEntries();
128 // Start loop on tracks in the hits containers
130 for (Int_t track=0; track<ntracks;track++) {
131 printf ("Processing Track: %d\n",track);
133 Int_t nbytes += TH->GetEvent(track);
135 TClonesArray *PadHits = RICH->PadHits(); // Cluster branch
136 TClonesArray *Hits = RICH->Hits(); // Hits branch
137 TClonesArray *Cerenkovs = RICH->Cerenkovs(); // Cerenkovs branch
139 //see hits distribution
140 Int_t nhits = Hits->GetEntriesFast();
141 if (nhits) Nh+=nhits;
142 for (Int_t hit=0;hit<nhits;hit++) {
143 mHit = (AliRICHHit*) Hits->UncheckedAt(hit);
144 Int_t nch = mHit->fChamber; // chamber number
145 Float_t x = mHit->fX; // x-pos of hit
146 Float_t y = mHit->fZ; // y-pos
147 Float_t z = mHit->fY;
148 Int_t index = mHit->fTrack;
149 Int_t particle = mHit->fParticle;
152 TParticle *current = (TParticle*)(*gAlice->Particles())[index];
154 Float_t energy=current->Energy();
156 R=TMath::Sqrt(current->Vx()*current->Vx() + current->Vy()*current->Vy());
158 //printf("Particle type: %d\n",current->GetPdgCode());
159 if (TMath::Abs(particle) < 50000051)
161 //if (TMath::Abs(particle) == 50000050 || TMath::Abs(particle) == 2112)
162 if (TMath::Abs(particle) == 2112 || TMath::Abs(particle) == 50000050)
164 //gMC->Rndm(&random, 1);
165 if (random->Rndm() < .1)
166 production->Fill(current->Vz(),R,(float) 1);
167 if (TMath::Abs(particle) == 50000050)
173 if (current->Energy()>0.001)
174 highprimaryphotons +=1;
177 if (TMath::Abs(particle) == 2112)
180 if (current->Energy()>0.0001)
186 production->Fill(current->Vz(),R,(float) 1);
187 //printf("Adding %d at %f\n",particle,R);
189 //mip->Fill(x,y,(float) 1);
192 if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
194 pionspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
195 if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
196 pionspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
199 pionspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
200 //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\R:%f\n\n\n\n\n\n\n\n\n",R);
202 //printf("Pion mass: %e\n",current->GetCalcMass());
204 if (TMath::Abs(particle)==211)
210 if (current->Energy()>1)
211 highprimarypions +=1;
215 if (TMath::Abs(particle)==2212)
217 protonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
218 if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
219 protonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
221 protonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
222 //printf("\n\n\n\n\n\n\nProton mass: %e\n\n\n\n\n\n\n\n\n",current->GetCalcMass());
225 if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
226 || TMath::Abs(particle)==311)
228 kaonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
229 if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
230 kaonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
232 kaonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
233 //printf("Kaon mass: %e\n",current->GetCalcMass());
235 if (TMath::Abs(particle)==321)
241 if (current->Energy()>1)
242 highprimarykaons +=1;
246 if (TMath::Abs(particle)==11)
248 electronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
249 if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
250 electronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
252 electronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
253 //printf("Electron mass: %e\n",current->GetCalcMass());
259 if (TMath::Abs(particle)==13)
261 muonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
262 if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
263 muonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
265 muonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
266 //printf("Muon mass: %e\n",current->GetCalcMass());
269 if (TMath::Abs(particle)==2112)
271 neutronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
272 if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
273 neutronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
276 neutronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
277 //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\R:%f\n\n\n\n\n\n\n\n\n",R);
279 //printf("Neutron mass: %e\n",current->GetCalcMass());
282 if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
284 if (current->Energy()-current->GetCalcMass()>1)
286 chargedspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
287 if (current->Vx()>.005 && current->Vy()>.005 && current->Vz()>.005)
288 chargedspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
290 chargedspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
293 //printf("Hits:%d\n",hit);
294 //printf ("Chamber number:%d x:%f y:%f\n",nch,x,y);
295 // Fill the histograms
297 //h->Fill(x,y,(float) 1);
306 //Create canvases, set the view range, show histograms
308 TCanvas *c15 = new TCanvas("c15","Mothers Production Vertices",50,50,600,600);
309 production->SetFillColor(42);
310 production->SetXTitle("z (m)");
311 production->SetYTitle("R (m)");
314 TCanvas *c16 = new TCanvas("c16","Particles Spectra II",100,100,600,350);
318 //TCanvas *c13 = new TCanvas("c13","Electron Spectra",400,10,600,700);
319 electronspectra1->SetFillColor(42);
320 electronspectra1->SetXTitle("log(GeV)");
321 electronspectra2->SetFillColor(46);
322 electronspectra2->SetXTitle("log(GeV)");
323 electronspectra3->SetFillColor(10);
324 electronspectra3->SetXTitle("log(GeV)");
326 electronspectra1->Draw();
327 electronspectra2->Draw("same");
328 electronspectra3->Draw("same");
331 //TCanvas *c14 = new TCanvas("c14","Muon Spectra",400,10,600,700);
332 muonspectra1->SetFillColor(42);
333 muonspectra1->SetXTitle("log(GeV)");
334 muonspectra2->SetFillColor(46);
335 muonspectra2->SetXTitle("log(GeV)");
336 muonspectra3->SetFillColor(10);
337 muonspectra3->SetXTitle("log(GeV)");
339 muonspectra1->Draw();
340 muonspectra2->Draw("same");
341 muonspectra3->Draw("same");
344 //TCanvas *c16 = new TCanvas("c16","Neutron Spectra",400,10,600,700);
345 //neutronspectra1->SetFillColor(42);
346 //neutronspectra1->SetXTitle("log(GeV)");
347 //neutronspectra2->SetFillColor(46);
348 //neutronspectra2->SetXTitle("log(GeV)");
349 //neutronspectra3->SetFillColor(10);
350 //neutronspectra3->SetXTitle("log(GeV)");
352 //neutronspectra1->Draw();
353 //neutronspectra2->Draw("same");
354 //neutronspectra3->Draw("same");
356 TCanvas *c9 = new TCanvas("c9","Particles Spectra",150,150,600,700);
357 //TCanvas *c9 = new TCanvas("c9","Pion Spectra",400,10,600,700);
361 pionspectra1->SetFillColor(42);
362 pionspectra1->SetXTitle("log(GeV)");
363 pionspectra2->SetFillColor(46);
364 pionspectra2->SetXTitle("log(GeV)");
365 pionspectra3->SetFillColor(10);
366 pionspectra3->SetXTitle("log(GeV)");
368 pionspectra1->Draw();
369 pionspectra2->Draw("same");
370 pionspectra3->Draw("same");
373 //TCanvas *c10 = new TCanvas("c10","Proton Spectra",400,10,600,700);
374 protonspectra1->SetFillColor(42);
375 protonspectra1->SetXTitle("log(GeV)");
376 protonspectra2->SetFillColor(46);
377 protonspectra2->SetXTitle("log(GeV)");
378 protonspectra3->SetFillColor(10);
379 protonspectra3->SetXTitle("log(GeV)");
381 protonspectra1->Draw();
382 protonspectra2->Draw("same");
383 protonspectra3->Draw("same");
386 //TCanvas *c11 = new TCanvas("c11","Kaon Spectra",400,10,600,700);
387 kaonspectra1->SetFillColor(42);
388 kaonspectra1->SetXTitle("log(GeV)");
389 kaonspectra2->SetFillColor(46);
390 kaonspectra2->SetXTitle("log(GeV)");
391 kaonspectra3->SetFillColor(10);
392 kaonspectra3->SetXTitle("log(GeV)");
394 kaonspectra1->Draw();
395 kaonspectra2->Draw("same");
396 kaonspectra3->Draw("same");
399 //TCanvas *c12 = new TCanvas("c12","Charged Particles Spectra",400,10,600,700);
400 chargedspectra1->SetFillColor(42);
401 chargedspectra1->SetXTitle("log(GeV)");
402 chargedspectra2->SetFillColor(46);
403 chargedspectra2->SetXTitle("log(GeV)");
404 chargedspectra3->SetFillColor(10);
405 chargedspectra3->SetXTitle("log(GeV)");
407 chargedspectra1->Draw();
408 chargedspectra2->Draw("same");
409 chargedspectra3->Draw("same");
411 // calculate the number of pads which give a signal
415 for (Int_t i=0; i< NpadX;i++) {
416 for (Int_t j=0;j< NpadY;j++) {
422 //printf("The total number of pads which give a signal: %d %d\n",Nh,Nh1);
424 printf("*****************************************\n");
425 printf("* Particle * Flux(m2) *\n");
426 printf("*****************************************\n");
428 printf("* Pions: * %3.1f *\n",pion/11.757312);
429 printf("* Charged Pions: * %3.1f *\n",chargedpions/11.757312);
430 printf("* Primary Pions: * %3.1f *\n",primarypions/11.757312);
431 printf("* Primary Pions (p>1GeV/c): * %3.1f *\n",highprimarypions/11.757312);
432 printf("* Kaons: * %3.1f *\n",kaon/11.757312);
433 printf("* Charged Kaons: * %3.1f *\n",chargedkaons/11.757312);
434 printf("* Primary Kaons: * %3.1f *\n",primarykaons/11.757312);
435 printf("* Primary Kaons (p>1GeV/c): * %3.1f *\n",highprimarykaons/11.757312);
436 printf("* Muons: * %3.1f *\n",muon/11.757312);
437 printf("* Electrons: * %3.1f *\n",electron/11.757312);
438 printf("* Positrons: * %3.1f *\n",positron/11.757312);
439 printf("* Protons: * %3.1f *\n",proton/11.757312);
440 printf("* All Charged: * %3.1f *\n",(chargedpions+chargedkaons+muon+electron+positron+proton)/11.757312);
441 printf("*****************************************\n");
442 //printf("* Photons: * %3.1f *\n",photons/11.757312);
443 //printf("* Primary Photons: * %3.1f *\n",primaryphotons/11.757312);
444 //printf("* Primary Photons (p>1MeV/c):* %3.1f *\n",highprimaryphotons/11.757312);
445 //printf("*****************************************\n");
446 //printf("* Neutrons: * %3.1f *\n",neutron);
447 //printf("* Neutrons (p>100keV/c): * %3.1f *\n",highneutrons);
448 //printf("*****************************************\n");
450 printf("\nEnd of macro\n");