]> git.uio.no Git - u/mrichter/AliRoot.git/blame - RICH/RICHpadtest.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Dummy copy constructor added
[u/mrichter/AliRoot.git] / RICH / RICHpadtest.C
CommitLineData
5bd23319 1void RICHpadtest (Int_t diaglevel,Int_t evNumber1=0,Int_t evNumber2=0)
2{
d2269328 3
5bd23319 4// Int_t diaglevel=3; // 1->Hits, 2->Spectra, 3->Statistics
d2269328 5
ddae0931 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
d2269328 14 Int_t NpadX = 162; // number of pads on X
15 Int_t NpadY = 162; // number of pads on Y
ddae0931 16
d2269328 17 Int_t Pad[162][162];
ddae0931 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 }
d2269328 23 gClassTable->GetID("AliRun");
ddae0931 24
25
26// Dynamically link some shared libs
27
28 if (gClassTable->GetID("AliRun") < 0) {
d2269328 29 gROOT->LoadMacro("loadlibs.C");
30 loadlibs();
31 }
32 else {
33 //delete gAlice;
34 gAlice = 0;
ddae0931 35 }
ddae0931 36
d2269328 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");
43
ddae0931 44// Get AliRun object from file or create it if not on file
d2269328 45
46 if (!gAlice) {
ddae0931 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");
d2269328 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 }
ddae0931 57
58// Create some histograms
59
d2269328 60 Int_t xmin= -NpadX/2;
61 Int_t xmax= NpadX/2;
62 Int_t ymin= -NpadY/2;
63 Int_t ymax= NpadY/2;
ddae0931 64
5bd23319 65 TH2F *hc0 = new TH2F("hc0","Zoom on center of Chamber 1",100,-50,50,100,-50,50);
d2269328 66 TH2F *hc1 = new TH2F("hc1","Chamber 1 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
67 TH2F *hc2 = new TH2F("hc2","Chamber 2 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
68 TH2F *hc3 = new TH2F("hc3","Chamber 3 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
69 TH2F *hc4 = new TH2F("hc4","Chamber 4 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
70 TH2F *hc5 = new TH2F("hc5","Chamber 5 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
71 TH2F *hc6 = new TH2F("hc6","Chamber 6 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
72 TH2F *hc7 = new TH2F("hc7","Chamber 7 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
73 TH2F *h = new TH2F("h","Detector hit distribution",150,-300,300,150,-300,300);
74 TH1F *Clcharge = new TH1F("Clcharge","Cluster Charge Distribution",500,0.,500.);
75 TH2F *cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-300,300,150,-300,300);
76 TH1F *ckovangle = new TH1F("ckovangle","Cerenkov angle per photon",200,.5,1);
77 TH1F *hckphi = new TH1F("hckphi","Cerenkov phi angle per photon",620,-3.1,3.1);
78 TH2F *feedback = new TH2F("feedback","Feedback hit distribution",150,-300,300,150,-300,300);
79 TH2F *mip = new TH2F("mip","Mip hit distribution",150,-300,300,150,-300,300);
80 TH1F *mother = new TH1F("mother","Cerenkovs per Mip",75,0.,75.);
81 TH1F *radius = new TH1F("radius","Mean distance to Mip",100,0.,20.);
82 TH1F *phspectra1 = new TH1F("phspectra","Photon Spectra",200,5.,10.);
83 TH1F *phspectra2 = new TH1F("phspectra","Photon Spectra",200,5.,10.);
84 TH1F *totalphotons = new TH1F("totalphotons","Produced Photons per Mip",100,200,700.);
85 TH1F *feedbacks = new TH1F("feedbacks","Produced Feedbacks per Mip",50,0.5,50.);
86 TH1F *padnumber = new TH1F("padnumber","Number of pads per cluster",50,-0.5,50.);
87 TH1F *padsev = new TH1F("padsev","Number of pads hit per event",50,0.5,100.);
88 TH1F *clusev = new TH1F("clusev","Number of clusters per event",50,0.5,50.);
89 TH1F *photev = new TH1F("photev","Number of photons per event",50,0.5,50.);
90 TH1F *feedev = new TH1F("feedev","Number of feedbacks per event",50,0.5,50.);
91 TH1F *padsmip = new TH1F("padsmip","Number of pads per event inside MIP region",50,0.5,50.);
92 TH1F *padscl = new TH1F("padscl","Number of pads per event from cluster count",50,0.5,100.);
93 TH1F *pionspectra = new TH1F("pionspectra","Pion Spectra",200,.5,10.);
94 TH1F *protonspectra = new TH1F("protonspectra","Proton Spectra",200,.5,10.);
95 TH1F *kaonspectra = new TH1F("kaonspectra","Kaon Spectra",100,.5,10.);
96 TH1F *kaonspectra = new TH1F("kaonspectra","Kaon Spectra",100,.5,10.);
97 TH1F *chargedspectra = new TH1F("chargedspectra","Charged particles above 1 GeV Spectra",100,.5,10.);
98 TH1F *hitsX = new TH1F("digitsX","Distribution of hits along x-axis",200,-300,300);
99 TH1F *hitsY = new TH1F("digitsY","Distribution of hits along z-axis",200,-300,300);
100 TH1F *hitsPhi = new TH1F("hitsPhi","Distribution of phi angle of incidence",100,-180,180);
101 TH1F *hitsTheta = new TH1F("hitsTheta","Distribution of Theta angle of incidence",100,0,15);
102 TH1F *Omega = new TH1F("omega","Reconstructed Cerenkov angle per track",200,.5,1);
103 TH1F *Theta = new TH1F("theta","Reconstructed theta incidence angle per track",200,0,15);
104 TH1F *Phi = new TH1F("phi","Reconstructed phi incidence per track",200,-180,180);
105
106
107
ddae0931 108
109// Start loop over events
110
111 Int_t Nh=0;
d2269328 112 Int_t pads=0;
ddae0931 113 Int_t Nh1=0;
d2269328 114 Int_t mothers[80000];
115 Int_t mothers2[80000];
116 Float_t mom[3];
117 Int_t nraw=0;
118 Int_t phot=0;
119 Int_t feed=0;
120 Int_t padmip=0;
121 for (Int_t i=0;i<100;i++) mothers[i]=0;
ddae0931 122 for (int nev=0; nev<= evNumber2; nev++) {
d2269328 123 Int_t nparticles = gAlice->GetEvent(nev);
124
125
126 //cout<<"nev "<<nev<<endl;
5bd23319 127 printf ("\nProcessing event: %d\n",nev);
d2269328 128 //cout<<"nparticles "<<nparticles<<endl;
129 printf ("Particles : %d\n",nparticles);
130 if (nev < evNumber1) continue;
131 if (nparticles <= 0) return;
132
ddae0931 133// Get pointers to RICH detector and Hits containers
d2269328 134
135 AliRICH *RICH = (AliRICH*)gAlice->GetDetector("RICH");
136 Int_t nent=(Int_t)gAlice->TreeR()->GetEntries();
137 gAlice->TreeR()->GetEvent(nent-1);
138 TClonesArray *Rawclusters = RICH->RawClustAddress(2); // Raw clusters branch
139 //printf ("Rawclusters:%p",Rawclusters);
140 Int_t nrawclusters = Rawclusters->GetEntriesFast();
141 //printf (" nrawclusters:%d\n",nrawclusters);
142 gAlice->TreeR()->GetEvent(nent-1);
143 TClonesArray *RecHits = RICH->RecHitsAddress(2);
144 Int_t nrechits = RecHits->GetEntriesFast();
145 //printf (" nrechits:%d\n",nrechits);
146 TTree *TH = gAlice->TreeH();
147 Int_t ntracks = TH->GetEntries();
ddae0931 148
ddae0931 149
d2269328 150
151 Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
152 gAlice->TreeD()->GetEvent(nent-1);
153
154
155 TClonesArray *Digits = RICH->DigitsAddress(2); // Raw clusters branch
156 Int_t ndigits = Digits->GetEntriesFast();
157 //printf("Digits:%d\n",ndigits);
158 padsev->Fill(ndigits,(float) 1);
159
160 for (Int_t ich=0;ich<7;ich++)
161 {
162 TClonesArray *Digits = RICH->DigitsAddress(ich); // Raw clusters branch
163 Int_t ndigits = Digits->GetEntriesFast();
164 //printf("Digits:%d\n",ndigits);
165 padsev->Fill(ndigits,(float) 1);
166 if (ndigits) {
167 for (Int_t hit=0;hit<ndigits;hit++) {
168 dHit = (AliRICHDigit*) Digits->UncheckedAt(hit);
169 //Int_t nchamber = padHit->fChamber; // chamber number
170 //Int_t nhit = dHit->fHitNumber; // hit number
171 Int_t qtot = dHit->fSignal; // charge
172 Int_t ipx = dHit->fPadX; // pad number on X
173 Int_t ipy = dHit->fPadY; // pad number on Y
174 //Int_t iqpad = dHit->fQpad; // charge per pad
175 //Int_t rpad = dHit->fRSec; // R-position of pad
5bd23319 176 //printf ("Pad hit, PadX:%d, PadY:%d\n",ipx,ipy);
177 if( ipx<=100 && ipy <=100 && ich==0) hc0->Fill(ipx,ipy,(float) qtot);
d2269328 178 if( ipx<=162 && ipy <=162 && ich==0) hc1->Fill(ipx,ipy,(float) qtot);
179 if( ipx<=162 && ipy <=162 && ich==1) hc2->Fill(ipx,ipy,(float) qtot);
180 if( ipx<=162 && ipy <=162 && ich==2) hc3->Fill(ipx,ipy,(float) qtot);
181 if( ipx<=162 && ipy <=162 && ich==3) hc4->Fill(ipx,ipy,(float) qtot);
182 if( ipx<=162 && ipy <=162 && ich==4) hc5->Fill(ipx,ipy,(float) qtot);
183 if( ipx<=162 && ipy <=162 && ich==5) hc6->Fill(ipx,ipy,(float) qtot);
184 if( ipx<=162 && ipy <=162 && ich==6) hc7->Fill(ipx,ipy,(float) qtot);
185 }
186 }
187 }
188
ddae0931 189// Start loop on tracks in the hits containers
d2269328 190 Int_t Nc=0;
191 for (Int_t track=0; track<ntracks;track++) {
192 printf ("Processing Track: %d\n",track);
193 gAlice->ResetHits();
194 Int_t nbytes += TH->GetEvent(track);
195 if (RICH) {
196 //RICH->ResetRawClusters();
197 TClonesArray *PadHits = RICH->PadHits(); // Cluster branch
198 TClonesArray *Hits = RICH->Hits(); // Hits branch
199 TClonesArray *Cerenkovs = RICH->Cerenkovs(); // Cerenkovs branch
200 }
201 //see hits distribution
202
ddae0931 203
d2269328 204 Int_t nhits = Hits->GetEntriesFast();
205 if (nhits) Nh+=nhits;
206 //printf("nhits %d\n",nhits);
207 for (Int_t hit=0;hit<nhits;hit++) {
208 mHit = (AliRICHHit*) Hits->UncheckedAt(hit);
ddae0931 209 Int_t nch = mHit->fChamber; // chamber number
210 Float_t x = mHit->fX; // x-pos of hit
d2269328 211 Float_t y = mHit->fZ; // y-pos
212 Float_t phi = mHit->fPhi; //Phi angle of incidence
213 Float_t theta = mHit->fTheta; //Theta angle of incidence
214 Int_t index = mHit->fTrack;
215 Int_t particle = mHit->fParticle;
216 Int_t freon = mHit->fLoss;
217
218 hitsX->Fill(x,(float) 1);
219 hitsY->Fill(y,(float) 1);
220
221 //printf("Particle:%d\n",particle);
222
223 TParticle *current = (TParticle*)(*gAlice->Particles())[index];
224 //printf("Particle type: %d\n",current->GetPdgCode());
225 if (TMath::Abs(particle) < 50000000)
226 {
227 mip->Fill(x,y,(float) 1);
228 if (current->Energy() - current->GetCalcMass()>1 && freon==1)
229 {
230 hitsPhi->Fill(phi,(float) 1);
231 hitsTheta->Fill(theta,(float) 1);
232 //printf("Theta:%f, Phi:%f\n",theta,phi);
233 }
234 }
235
236 if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
237 {
238 pionspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
239 }
240 if (TMath::Abs(particle)==2212)
241 {
242 protonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
243 }
244 if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
245 || TMath::Abs(particle)==311)
246 {
247 kaonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
248 }
249 if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
250 {
251 if (current->Energy() - current->GetCalcMass()>1)
252 chargedspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
253 }
254 //printf("Hits:%d\n",hit);
255 //printf ("Chamber number:%d x:%f y:%f\n",nch,x,y);
ddae0931 256 // Fill the histograms
d2269328 257 Nh1+=nhits;
258 h->Fill(x,y,(float) 1);
259 //}
260 //}
ddae0931 261 }
d2269328 262
263 Int_t ncerenkovs = Cerenkovs->GetEntriesFast();
264
265 if (ncerenkovs) {
266 for (Int_t hit=0;hit<ncerenkovs;hit++) {
267 cHit = (AliRICHCerenkov*) Cerenkovs->UncheckedAt(hit);
268 Int_t nchamber = cHit->fChamber; // chamber number
269 Int_t index = cHit->fTrack;
270 Int_t pindex = cHit->fIndex;
271 Int_t cx = cHit->fX; // x-position
272 Int_t cy = cHit->fZ; // y-position
273 Int_t cmother = cHit->fCMother; // Index of mother particle
274 Int_t closs = cHit->fLoss; // How did the particle get lost?
275 //printf ("Cerenkov hit, X:%d, Y:%d\n",cx,cy);
276
277
278
279 TParticle *current = (TParticle*)(*gAlice->Particles())[index];
280
281 if (current->GetPdgCode() == 50000051)
282 {
283 if (closs==4)
284 {
285 feedback->Fill(cx,cy,(float) 1);
286 feed++;
287 }
288 }
289 if (current->GetPdgCode() == 50000050)
290 {
291 if (closs==4)
292 {
293 cerenkov->Fill(cx,cy,(float) 1);
294
295
296
297 TParticle *MIP = (TParticle*)(*gAlice->Particles())[cmother];
298 mipHit = (AliRICHHit*) Hits->UncheckedAt(0);
299 mom[0] = current->Px();
300 mom[1] = current->Py();
301 mom[2] = current->Pz();
302 Float_t energyckov = current->Energy();
303 /*mom[0] = cHit->fMomX;
304 mom[1] = cHit->fMomZ;
305 mom[2] = cHit->fMomY;*/
306 Float_t energymip = MIP->Energy();
307 Float_t Mip_px = mipHit->fMomX;
308 Float_t Mip_py = mipHit->fMomY;
309 Float_t Mip_pz = mipHit->fMomZ;
310 /*Float_t Mip_px = MIP->Px();
311 Float_t Mip_py = MIP->Py();
312 Float_t Mip_pz = MIP->Pz();*/
313
314
315
316 Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
317 Float_t rt = TMath::Sqrt(r);
318 Float_t Mip_r = Mip_px*Mip_px + Mip_py*Mip_py + Mip_pz*Mip_pz;
319 Float_t Mip_rt = TMath::Sqrt(Mip_r);
320 Float_t coscerenkov = (mom[0]*Mip_px + mom[1]*Mip_py + mom[2]*Mip_pz)/(rt*Mip_rt);
321 Float_t cherenkov = TMath::ACos(coscerenkov);
322 ckovangle->Fill(cherenkov,(float) 1); //Cerenkov angle calculus
323 //printf("Cherenkov: %f\n",cherenkov);
324 Float_t ckphi=TMath::ATan2(mom[0], mom[2]);
325 hckphi->Fill(ckphi,(float) 1);
326
327
328 Float_t mix = MIP->Vx();
329 Float_t miy = MIP->Vy();
330 Float_t mx = mipHit->fX;
331 Float_t my = mipHit->fZ;
332 //printf("FX %e, FY %e, VX %e, VY %e\n",cx,cy,mx,my);
333 Float_t dx = cx - mx;
334 Float_t dy = cy - my;
335 //printf("Dx:%f, Dy:%f\n",dx,dy);
336 Float_t final_radius = TMath::Sqrt(dx*dx+dy*dy);
337 //printf("Final radius:%f\n",final_radius);
338 radius->Fill(final_radius,(float) 1);
339
340 phspectra1->Fill(energyckov*1e9,(float) 1);
341 phot++;
342 }
343 else
344 phspectra2->Fill(energyckov*1e9,(float) 1);
345 for (Int_t nmothers=0;nmothers<=ntracks;nmothers++){
346 if (cmother == nmothers){
347 if (closs == 4)
348 mothers2[cmother]++;
349 mothers[cmother]++;
350 }
351 }
352 }
353 }
354 }
355
356 if (nrawclusters) {
357 for (Int_t hit=0;hit<nrawclusters;hit++) {
358 rcHit = (AliRICHRawCluster*) Rawclusters->UncheckedAt(hit);
359 //Int_t nchamber = rcHit->fChamber; // chamber number
360 //Int_t nhit = cHit->fHitNumber; // hit number
361 Int_t qtot = rcHit->fQ; // charge
362 Int_t fx = rcHit->fX; // x-position
363 Int_t fy = rcHit->fY; // y-position
364 Int_t type = rcHit->fCtype; // cluster type ?
365 Int_t mult = rcHit->fMultiplicity; // How many pads form the cluster
366 pads += mult;
367 if (qtot > 0) {
368 if (fx>-4 && fx<4 && fy>-4 && fy<4) {
369 padmip+=mult;
370 } else {
371 padnumber->Fill(mult,(float) 1);
372 nraw++;
373 if (mult<4) Clcharge->Fill(qtot,(float) 1);
374 }
375 }
376 }
377 }
378
379 if(nrechits)
380 {
381 for (Int_t hit=0;hit<nrechits;hit++) {
382 recHit = (AliRICHRecHit*) RecHits->UncheckedAt(hit);
383 Float_t r_omega = recHit->Omega; // Cerenkov angle
384 Float_t r_theta = recHit->Theta; // Theta angle of incidence
385 Float_t r_phi = recHit->Phi; // Phi angle if incidence
386
387 Omega->Fill(r_omega,(float) 1);
388 Theta->Fill(r_theta*180/TMath::Pi(),(float) 1);
389 Phi->Fill(r_phi*180/TMath::Pi(),(float) 1);
390
391 //printf("Omega: %f, Theta: %f, Phi: %f\n",r_omega,r_theta,r_phi);
392 }
393 }
ddae0931 394 }
d2269328 395
396 for (Int_t nmothers=0;nmothers<ntracks;nmothers++){
397 totalphotons->Fill(mothers[nmothers],(float) 1);
398 mother->Fill(mothers2[nmothers],(float) 1);
399 //printf ("Entries in %d : %d\n",nmothers, mothers[nmothers]);
400 }
401
402 clusev->Fill(nraw,(float) 1);
403 photev->Fill(phot,(float) 1);
404 feedev->Fill(feed,(float) 1);
405 padsmip->Fill(padmip,(float) 1);
406 padscl->Fill(pads,(float) 1);
407 //printf("Photons:%d\n",phot);
408 phot = 0;
409 feed = 0;
410 pads = 0;
411 nraw=0;
412 padmip=0;
ddae0931 413
d2269328 414 }
ddae0931 415
d2269328 416 //Create canvases, set the view range, show histograms
ddae0931 417
d2269328 418 switch(diaglevel)
419 {
420 case 1:
ddae0931 421
d2269328 422 if (ndigits)
423 {
424 TCanvas *c1 = new TCanvas("c1","Alice RICH pad hits",50,10,1200,700);
425 c1->Divide(4,2);
426 c1->cd(1);
427 hc1->SetXTitle("ix (npads)");
428 hc1->Draw();
429 c1->cd(2);
430 hc2->SetXTitle("ix (npads)");
431 hc2->Draw();
432 c1->cd(3);
433 hc3->SetXTitle("ix (npads)");
434 hc3->Draw();
435 c1->cd(4);
436 hc4->SetXTitle("ix (npads)");
437 hc4->Draw();
438 c1->cd(5);
439 hc5->SetXTitle("ix (npads)");
440 hc5->Draw();
441 c1->cd(6);
442 hc6->SetXTitle("ix (npads)");
443 hc6->Draw();
444 c1->cd(7);
445 hc7->SetXTitle("ix (npads)");
446 hc7->Draw();
5bd23319 447 c1->cd(8);
448 hc0->SetXTitle("ix (npads)");
449 hc0->Draw();
d2269328 450 }
451//
452 TCanvas *c4 = new TCanvas("c4","Hits per type",400,10,600,700);
453 c4->Divide(2,2);
454
455 c4->cd(1);
456 feedback->SetFillColor(42);
457 feedback->SetXTitle("x (pads)");
458 feedback->SetYTitle("y (pads)");
459 feedback->Draw();
460
461 c4->cd(2);
462 mip->SetFillColor(42);
463 mip->SetXTitle("x (pads)");
464 mip->SetYTitle("y (pads)");
465 mip->Draw();
466
467 c4->cd(3);
468 cerenkov->SetFillColor(42);
469 cerenkov->SetXTitle("x (pads)");
470 cerenkov->SetYTitle("y (pads)");
471 cerenkov->Draw();
472
473 c4->cd(4);
474 h->SetFillColor(42);
475 h->SetXTitle("x (cm)");
476 h->SetYTitle("y (cm)");
477 h->Draw();
478
479 TCanvas *c10 = new TCanvas("c10","Hits distribution",400,10,600,350);
480 c10->Divide(2,1);
481
482 c10->cd(1);
483 hitsX->SetFillColor(42);
5bd23319 484 hitsX->SetXTitle("(cm)");
d2269328 485 hitsX->Draw();
486
487 c10->cd(2);
488 hitsY->SetFillColor(42);
5bd23319 489 hitsY->SetXTitle("(cm)");
d2269328 490 hitsY->Draw();
491
492
493 break;
494//
495 case 2:
496
497 TCanvas *c6 = new TCanvas("c6","Photon Spectra",50,10,600,350);
498 c6->Divide(2,1);
499
500 c6->cd(1);
501 phspectra2->SetFillColor(42);
502 phspectra2->SetXTitle("energy (eV)");
503 phspectra2->Draw();
504 c6->cd(2);
505 phspectra1->SetFillColor(42);
506 phspectra1->SetXTitle("energy (eV)");
507 phspectra1->Draw();
508
509 TCanvas *c9 = new TCanvas("c9","Particles Spectra",400,10,600,700);
510 c9->Divide(2,2);
511
512 c9->cd(1);
513 pionspectra->SetFillColor(42);
514 pionspectra->SetXTitle("(GeV)");
515 pionspectra->Draw();
516
517 c9->cd(2);
518 protonspectra->SetFillColor(42);
519 protonspectra->SetXTitle("(GeV)");
520 protonspectra->Draw();
521
522 c9->cd(3);
523 kaonspectra->SetFillColor(42);
524 kaonspectra->SetXTitle("(GeV)");
525 kaonspectra->Draw();
526
527 c9->cd(4);
528 chargedspectra->SetFillColor(42);
529 chargedspectra->SetXTitle("(GeV)");
530 chargedspectra->Draw();
ddae0931 531
d2269328 532 break;
533
534 case 3:
535
536 if (nrawclusters) {
537 TCanvas *c3=new TCanvas("c3","Clusters Statistics",400,10,600,700);
538 c3->Divide(2,2);
539
540 c3->cd(1);
541 Clcharge->SetFillColor(42);
542 Clcharge->SetXTitle("ADC units");
543 Clcharge->Draw();
544
545 c3->cd(2);
546 padnumber->SetFillColor(42);
547 padnumber->SetXTitle("(counts)");
548 padnumber->Draw();
549
550 c3->cd(3);
551 clusev->SetFillColor(42);
552 clusev->SetXTitle("(counts)");
553 clusev->Draw();
ddae0931 554
d2269328 555 c3->cd(4);
556 padsmip->SetFillColor(42);
557 padsmip->SetXTitle("(counts)");
558 padsmip->Draw();
559 }
560
561 if (nev<1)
562 {
563 TCanvas *c11 = new TCanvas("c11","Cherenkov per Mip",400,10,600,700);
564 mother->SetFillColor(42);
565 mother->SetXTitle("counts");
566 mother->Draw();
567 }
568
569 TCanvas *c2 = new TCanvas("c2","Angles of incidence",50,10,600,700);
570 c2->Divide(2,2);
571
572 c2->cd(1);
573 hitsPhi->SetFillColor(42);
574 hitsPhi->Draw();
575 c2->cd(2);
576 hitsTheta->SetFillColor(42);
577 hitsTheta->Draw();
578 c2->cd(3);
579 Phi->SetFillColor(42);
580 Phi->Draw();
581 c2->cd(4);
582 Theta->SetFillColor(42);
583 Theta->Draw();
584
585
586 TCanvas *c5 = new TCanvas("c5","Ring Statistics",50,10,600,700);
587 c5->Divide(2,2);
588
589 c5->cd(1);
590 ckovangle->SetFillColor(42);
591 ckovangle->SetXTitle("angle (radians)");
592 ckovangle->Draw();
593
594 c5->cd(2);
595 radius->SetFillColor(42);
596 radius->SetXTitle("radius (cm)");
597 radius->Draw();
598
599 c5->cd(3);
600 Omega->SetFillColor(42);
601 Omega->SetXTitle("angle (radians)");
602 Omega->Draw();
603
604 TCanvas *c7 = new TCanvas("c7","Production Statistics",400,10,600,700);
605 c7->Divide(2,2);
606
607 c7->cd(1);
608 totalphotons->SetFillColor(42);
609 totalphotons->SetXTitle("Photons (counts)");
610 totalphotons->Draw();
611
612 c7->cd(2);
613 photev->SetFillColor(42);
614 photev->SetXTitle("(counts)");
615 photev->Draw();
616
617 c7->cd(3);
618 feedev->SetFillColor(42);
619 feedev->SetXTitle("(counts)");
620 feedev->Draw();
621
622 c7->cd(4);
623 padsev->SetFillColor(42);
624 padsev->SetXTitle("(counts)");
625 padsev->Draw();
626
627 break;
628
629 }
630
631
632 // calculate the number of pads which give a signal
633
634
635 Int_t Np=0;
636 for (Int_t i=0;i< NpadX;i++) {
637 for (Int_t j=0;j< NpadY;j++) {
638 if (Pad[i][j]>=6){
639 Np+=1;
640 }
641 }
642 }
643 //printf("The total number of pads which give a signal: %d %d\n",Nh,Nh1);
644 printf("End of macro\n");
ddae0931 645}
d2269328 646
647