]> git.uio.no Git - u/mrichter/AliRoot.git/blame - RICH/RICHpadtestC.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
add NaN trap to linux version
[u/mrichter/AliRoot.git] / RICH / RICHpadtestC.C
CommitLineData
6a748468 1Int_t diaglevel=2; // 1->Hits, 2->Spectra, 3->Statistics
2
3
4void RICHpadtestC (Int_t evNumber1=0,Int_t evNumber2=0)
5{
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.
10//
11/////////////////////////////////////////////////////////////////////////
12
13
14 Int_t NpadX = 162; // number of pads on X
15 Int_t NpadY = 162; // number of pads on Y
16
17 Int_t Pad[162][162];
18 for (Int_t i=0;i<NpadX;i++) {
19 for (Int_t j=0;j<NpadY;j++) {
20 Pad[i][j]=0;
21 }
22 }
23
24
25
26// Dynamically link some shared libs
27
28 if (gClassTable->GetID("AliRun") < 0) {
85f1fa61 29 gROOT->LoadMacro("loadlibs.C");
30 loadlibs();
31 }
32 else {
33 //delete gAlice;
34 gAlice = 0;
6a748468 35 }
6a748468 36
85f1fa61 37 gAlice=0;
38
39// Connect the Root Galice file containing Geometry, Kine and Hits
40
41 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
42 if (!file) file = new TFile("galice.root","UPDATE");
43
6a748468 44// Get AliRun object from file or create it if not on file
85f1fa61 45
46 if (!gAlice) {
6a748468 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");
85f1fa61 50 }
51 else {
52 delete gAlice;
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");
56 }
6a748468 57
58// Create some histograms
59
60 Int_t xmin= -NpadX/2;
61 Int_t xmax= NpadX/2;
62 Int_t ymin= -NpadY/2;
63 Int_t ymax= NpadY/2;
64
6a748468 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);
6a748468 86 TH2F *production = new TH2F("production","Mother production vertices",100,-300,300,100,0,600);
87
88
89
90
91// Start loop over events
92
93 Int_t Nh=0;
94 Int_t pads=0;
95 Int_t Nh1=0;
6a748468 96 Float_t mom[3];
6a748468 97 Int_t nraw=0;
98 Int_t phot=0;
99 Int_t feed=0;
100 Int_t padmip=0;
59f5e5cd 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;
104
105 TRandom random;
106
6a748468 107 for (int nev=0; nev<= evNumber2; nev++) {
108 Int_t nparticles = gAlice->GetEvent(nev);
109
110
6a748468 111 printf ("Event number : %d\n",nev);
6a748468 112 printf ("Number of particles: %d\n",nparticles);
113 if (nev < evNumber1) continue;
114 if (nparticles <= 0) return;
115
116// Get pointers to RICH detector and Hits containers
117
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
6a748468 122 Int_t nrawclusters = Rawclusters->GetEntriesFast();
6a748468 123 TTree *TH = gAlice->TreeH();
124 Int_t ntracks = TH->GetEntries();
125
126
127
6a748468 128// Start loop on tracks in the hits containers
129 Int_t Nc=0;
130 for (Int_t track=0; track<ntracks;track++) {
131 printf ("Processing Track: %d\n",track);
132 gAlice->ResetHits();
133 Int_t nbytes += TH->GetEvent(track);
134 if (RICH) {
6a748468 135 TClonesArray *PadHits = RICH->PadHits(); // Cluster branch
136 TClonesArray *Hits = RICH->Hits(); // Hits branch
137 TClonesArray *Cerenkovs = RICH->Cerenkovs(); // Cerenkovs branch
138 }
139 //see hits distribution
140 Int_t nhits = Hits->GetEntriesFast();
141 if (nhits) Nh+=nhits;
6a748468 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;
150 Float_t R;
151
6a748468 152 TParticle *current = (TParticle*)(*gAlice->Particles())[index];
153
85f1fa61 154 Float_t energy=current->Energy();
155
6a748468 156 R=TMath::Sqrt(current->Vx()*current->Vx() + current->Vy()*current->Vy());
157
158 //printf("Particle type: %d\n",current->GetPdgCode());
159 if (TMath::Abs(particle) < 50000051)
160 {
59f5e5cd 161 //if (TMath::Abs(particle) == 50000050 || TMath::Abs(particle) == 2112)
162 if (TMath::Abs(particle) == 2112 || TMath::Abs(particle) == 50000050)
6a748468 163 {
59f5e5cd 164 //gMC->Rndm(&random, 1);
165 if (random->Rndm() < .1)
6a748468 166 production->Fill(current->Vz(),R,(float) 1);
59f5e5cd 167 if (TMath::Abs(particle) == 50000050)
168 {
169 photons +=1;
170 if (R<.005)
171 {
172 primaryphotons +=1;
173 if (current->Energy()>0.001)
174 highprimaryphotons +=1;
175 }
176 }
177 if (TMath::Abs(particle) == 2112)
178 {
179 neutron +=1;
180 if (current->Energy()>0.0001)
181 highneutrons +=1;
182 }
6a748468 183 }
184 else
185 {
186 production->Fill(current->Vz(),R,(float) 1);
85f1fa61 187 //printf("Adding %d at %f\n",particle,R);
6a748468 188 }
189 //mip->Fill(x,y,(float) 1);
190 }
191
192 if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
193 {
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);
197 if (R>2.5 && R<4.5)
198 {
199 pionspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
85f1fa61 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);
6a748468 201 }
85f1fa61 202 //printf("Pion mass: %e\n",current->GetCalcMass());
6a748468 203 pion +=1;
59f5e5cd 204 if (TMath::Abs(particle)==211)
205 {
206 chargedpions +=1;
207 if (R<.005)
208 {
209 primarypions +=1;
210 if (current->Energy()>1)
211 highprimarypions +=1;
212 }
213 }
6a748468 214 }
215 if (TMath::Abs(particle)==2212)
216 {
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);
220 if (R>3 && R<4.3)
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());
223 proton +=1;
224 }
225 if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
226 || TMath::Abs(particle)==311)
227 {
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);
231 if (R>2.5 && R<4.5)
232 kaonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
85f1fa61 233 //printf("Kaon mass: %e\n",current->GetCalcMass());
6a748468 234 kaon +=1;
59f5e5cd 235 if (TMath::Abs(particle)==321)
236 {
237 chargedkaons +=1;
238 if (R<.005)
239 {
240 primarykaons +=1;
241 if (current->Energy()>1)
242 highprimarykaons +=1;
243 }
244 }
6a748468 245 }
246 if (TMath::Abs(particle)==11)
247 {
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);
251 if (R>2.5 && R<4.5)
252 electronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
85f1fa61 253 //printf("Electron mass: %e\n",current->GetCalcMass());
59f5e5cd 254 if (particle == 11)
85f1fa61 255 electron +=1;
256 if (particle == -11)
59f5e5cd 257 positron +=1;
6a748468 258 }
259 if (TMath::Abs(particle)==13)
260 {
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);
264 if (R>2.5 && R<4.5)
265 muonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
85f1fa61 266 //printf("Muon mass: %e\n",current->GetCalcMass());
6a748468 267 muon +=1;
268 }
269 if (TMath::Abs(particle)==2112)
270 {
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);
274 if (R>2.5 && R<4.5)
275 {
276 neutronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
85f1fa61 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);
6a748468 278 }
85f1fa61 279 //printf("Neutron mass: %e\n",current->GetCalcMass());
6a748468 280 neutron +=1;
281 }
282 if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
283 {
284 if (current->Energy()-current->GetCalcMass()>1)
285 {
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);
289 if (R>2.5 && R<4.5)
290 chargedspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
291 }
292 }
293 //printf("Hits:%d\n",hit);
294 //printf ("Chamber number:%d x:%f y:%f\n",nch,x,y);
295 // Fill the histograms
296 Nh1+=nhits;
297 //h->Fill(x,y,(float) 1);
298 //}
299 //}
300 }
301
6a748468 302 }
303
6a748468 304 }
305
306 //Create canvases, set the view range, show histograms
307
85f1fa61 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)");
312 production->Draw();
6a748468 313
85f1fa61 314 TCanvas *c16 = new TCanvas("c16","Particles Spectra II",100,100,600,350);
315 c16->Divide(2,1);
316
317 c16->cd(1);
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)");
325 //c13->SetLogx();
326 electronspectra1->Draw();
327 electronspectra2->Draw("same");
328 electronspectra3->Draw("same");
329
330 c16->cd(2);
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)");
338 //c14->SetLogx();
339 muonspectra1->Draw();
340 muonspectra2->Draw("same");
341 muonspectra3->Draw("same");
342
343 //c16->cd(3);
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)");
351 //c16->SetLogx();
352 //neutronspectra1->Draw();
353 //neutronspectra2->Draw("same");
354 //neutronspectra3->Draw("same");
355
356 TCanvas *c9 = new TCanvas("c9","Particles Spectra",150,150,600,700);
357 //TCanvas *c9 = new TCanvas("c9","Pion Spectra",400,10,600,700);
358 c9->Divide(2,2);
359
360 c9->cd(1);
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)");
367 //c9->SetLogx();
368 pionspectra1->Draw();
369 pionspectra2->Draw("same");
370 pionspectra3->Draw("same");
371
372 c9->cd(2);
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)");
380 //c10->SetLogx();
381 protonspectra1->Draw();
382 protonspectra2->Draw("same");
383 protonspectra3->Draw("same");
384
385 c9->cd(3);
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)");
393 //c11->SetLogx();
394 kaonspectra1->Draw();
395 kaonspectra2->Draw("same");
396 kaonspectra3->Draw("same");
397
398 c9->cd(4);
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)");
406 //c12->SetLogx();
407 chargedspectra1->Draw();
408 chargedspectra2->Draw("same");
409 chargedspectra3->Draw("same");
410
6a748468 411 // calculate the number of pads which give a signal
412
413
414 Int_t Np=0;
415 for (Int_t i=0; i< NpadX;i++) {
416 for (Int_t j=0;j< NpadY;j++) {
417 if (Pad[i][j]>=6){
418 Np+=1;
419 }
420 }
421 }
422 //printf("The total number of pads which give a signal: %d %d\n",Nh,Nh1);
423
85f1fa61 424 printf("*****************************************\n");
425 printf("* Particle * Flux(m2) *\n");
426 printf("*****************************************\n");
427
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");
449
450 printf("\nEnd of macro\n");
6a748468 451}
452
453
454
455
59f5e5cd 456