Macro to find raw clusters (JB, AM)
[u/mrichter/AliRoot.git] / RICH / RICHtest.C
CommitLineData
ddae0931 1void RICHtest (Int_t evNumber1=0,Int_t evNumber2=0)
2{
3/////////////////////////////////////////////////////////////////////////
4// This macro is a small example of a ROOT macro
5// illustrating how to read the output of GALICE
6// and do some analysis.
7//
8/////////////////////////////////////////////////////////////////////////
9
10// Dynamically link some shared libs
11
12 if (gClassTable->GetID("AliRun") < 0) {
13 gROOT->LoadMacro("loadlibs.C");
14 loadlibs();
15 }
16// Connect the Root Galice file containing Geometry, Kine and Hits
17
18 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
19
20
21 if (!file) {
22 printf("\n Creating galice.root \n");
23 file = new TFile("galice.root");
24 } else {
25 printf("\n galice.root found in file list");
26 }
27 file->ls();
28
29// Get AliRun object from file or create it if not on file
30 if (!gAlice) {
31 gAlice = (AliRun*)(file->Get("gAlice"));
32 if (gAlice) printf("AliRun object found on file\n");
33 if (!gAlice) {
34 printf("\n create new gAlice object");
35 gAlice = new AliRun("gAlice","Alice test program");
36 }
37 }
38
39// Create some histograms
40
41
42 TH1F *ccharge1 = new TH1F("ccharge1","Cluster Charge distribution"
43 ,100,0.,500.);
44 TH1F *pcharge1 = new TH1F("pcharge1","Pad Charge distribution"
45 ,100,0.,200.);
46 TH1F *xresid1 = new TH1F("xresid1","x-Residuals"
47 ,100,-0.5,0.5);
48 TH1F *yresid1 = new TH1F("yresid1","y-Residuals"
49 ,100,-0.2,0.2);
50 TH1F *npads1 = new TH1F("npads1" ,"Pads per Hit"
51 ,20,-0.5,19.5);
52 TH2F *xresys1 = new TH2F("xresys1","x-Residuals systematics"
53 ,50,-2.0,2.0,100,-1.0,1.0);
54 TH2F *yresys1 = new TH2F("yresys1","y-Residuals systematics"
55 ,50,-0.4,0.4,100,-0.1,0.1);
56
57 TH1F *ccharge2 = new TH1F("ccharge2","Cluster Charge distribution"
58 ,100,0.,500.);
59 TH1F *pcharge2 = new TH1F("pcharge2","Pad Charge distribution"
60 ,100,0.,200.);
61 TH1F *xresid2 = new TH1F("xresid2","x-Residuals"
62 ,100,-0.5,0.5);
63 TH1F *yresid2 = new TH1F("yresid2","y-Residuals"
64 ,100,-0.2,0.2);
65 TH1F *npads2 = new TH1F("npads2" ,"Pads per Hit"
66 ,20,-0.5,19.5);
67 TH2F *xresys2 = new TH2F("xresys2","x-Residuals systematics"
68 ,50,-2.0,2.0,100,-1.0,1.0);
69 TH2F *yresys2 = new TH2F("yresys2","y-Residuals systematics"
70 ,50,-0.4,0.4,100,-0.1,0.1);
71
72
73 AliRICHchamber* iChamber;
74//
75// Loop over events
76//
77 Int_t Nh=0;
78 Int_t Nh1=0;
79 for (Int_t nev=0; nev<= evNumber2; nev++) {
80 //cout << "nev " << nev <<endl;
81 printf ("Event Number : %d\n",nev);
82 Int_t nparticles = gAlice->GetEvent(nev);
83 //cout << "nparticles " << nparticles <<endl;
84 printf ("Number of particles: %d\n",nparticles);
85 if (nev < evNumber1) continue;
86 if (nparticles <= 0) return;
87
88 AliRICH *RICH = (AliRICH*) gAlice->GetDetector("RICH");
89 printf("\n track %d %d \n ", nev, RICH);
90
91 TTree *TH = gAlice->TreeH();
92 Int_t ntracks = TH->GetEntries();
93 Int_t Nc=0;
94 Int_t npad[2];
95 Float_t Q[2], xmean[2],ymean[2],xres[2],yres[2], xonpad[2], yonpad[2];
96
97 printf ("Just got it");
98//
99// Loop over events
100//
101 for (Int_t track=0; track<ntracks;track++) {
102
103 gAlice->ResetHits();
104 Int_t nbytes += TH->GetEvent(track);
105 if (RICH) {
106 for(AliRICHhit* mHit=(AliRICHhit*)RICH->FirstHit(-1);
107 mHit;
108 mHit=(AliRICHhit*)RICH->NextHit())
109 {
110 Int_t nch = mHit->fChamber; // chamber number
111 Float_t x = mHit->fX; // x-pos of hit
112 Float_t y = mHit->fY; // y-pos
113//
114//
115 iChamber = & RICH->Chamber(nch-1);
116 response=iChamber->GetResponseModel(mip);
117//
118//
119 if (nch <= 1) {
120 for (Int_t i=0; i<2; i++) {
121 xmean[i]=0;
122 ymean[i]=0;
123 Q[i]=0;
124 npad[i]=0;
125 }
126 for (AliRICHcluster *mPad=(AliRICHcluster*)RICH
127 ->FirstPad(mHit);
128 mPad;
129 mPad=(AliRICHcluster*)RICH->NextPad()
130 )
131 {
132 Int_t nseg = mPad->fCathode; // segmentation
133 Int_t nhit = mPad->fHitNumber; // hit number
134 Int_t qtot = mPad->fQ; // charge
135 Int_t ipx = mPad->fPadX; // pad number on X
136 Int_t ipy = mPad->fPadY; // pad number on Y
137 Int_t iqpad = mPad->fQpad; // charge per pad
138 Int_t secpad = mPad->fRSec; // r-pos of pad
139//
140//
141 segmentation=iChamber->GetSegmentationModel(nseg);
142 Int_t ind=nseg-1;
143// printf("Pad of hit, padx, pady, iqpad, ncha %d %d %d %d %d %d\n",
144// nhit, ipx, ipy, iqpad, nseg, ind);
145
146
147 if (nseg == 1) {
148 pcharge1->Fill(iqpad,(float) 1);
149 ccharge1->Fill(qtot ,(float) 1);
150 } else {
151 pcharge2->Fill(iqpad,(float) 1);
152 ccharge2->Fill(qtot ,(float) 1);
153 }
154// Calculate centre of gravity
155//
156 if (iqpad == 0 && ind==0) {
157 printf("\n attention iqpad 0");
158 }
159
160 if (iqpad > 0 && secpad==1) {
161 Float_t xc,yc;
162 npad[ind]++;
163 segmentation->GetPadCxy(ipx,ipy,xc,yc);
164// printf("\n pad %d %d ", ipx, ipy);
165// printf("\n pos %f %f ", xc, yc);
166 xmean[ind]+=Float_t(iqpad*xc);
167 ymean[ind]+=Float_t(iqpad*yc);
168 Q[ind]+=iqpad;
169 }
170
171 } //pad hit loop
172 for (Int_t i=0; i<2; i++) {
173 segmentation = iChamber->GetSegmentationModel(i+1);
174
175// if (npad[i] ==0) {
176// printf("\n #pads=0 %f %f",x,y);
177// }
178
179 if (Q[i] >0) {
180 xmean[i] = xmean[i]/Q[i];
181 xres[i] = xmean[i]-x;
182 ymean[i] = ymean[i]/Q[i];
183 yres[i] = ymean[i]-y;
184// printf("\n xmean, ymean %f %f",xmean[i],ymean[i]);
185// printf("\n xhit, yhit %f %f",x,y);
186
187// Systematic Error
188//
189 Int_t icx, icy;
190 segmentation->
191 GetPadIxy(x,y,icx,icy);
192 segmentation->
193 GetPadCxy(icx,icy,xonpad[i],yonpad[i]);
194 xonpad[i]-=x;
195 yonpad[i]-=y;
196 }
197 } // plane loop
198 xresid1->Fill(xres[0],(float) 1);
199 yresid1->Fill(yres[0],(float) 1);
200 npads1->Fill(npad[0],(float) 1);
201 if (npad[0] >=2 && Q[0] > 6 ) {
202 xresys1->Fill(xonpad[0],xres[0],(float) 1);
203 yresys1->Fill(yonpad[0],yres[0],(float) 1);
204 }
205
206 xresid2->Fill(xres[1],(float) 1);
207 yresid2->Fill(yres[1],(float) 1);
208 npads2->Fill(npad[1],(float) 1);
209 if (npad[1] >=2 && Q[1] > 6) {
210 xresys2->Fill(xonpad[1],xres[1],(float) 1);
211 yresys2->Fill(yonpad[1],yres[1],(float) 1);
212 }
213 } // chamber 1
214 } // hit loop
215 } // if RICH
216 } // track loop
217 } // event loop
218
219//Create a canvas, set the view range, show histograms
220 TCanvas *c1 = new TCanvas("c1","Charge and Residuals",400,10,600,700);
221 pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
222 pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
223 pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
224 pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
225 pad11->SetFillColor(11);
226 pad12->SetFillColor(11);
227 pad13->SetFillColor(11);
228 pad14->SetFillColor(11);
229 pad11->Draw();
230 pad12->Draw();
231 pad13->Draw();
232 pad14->Draw();
233
234 pad11->cd();
235 ccharge1->SetFillColor(42);
236 ccharge1->SetXTitle("ADC units");
237 ccharge1->Draw();
238
239 pad12->cd();
240 pcharge1->SetFillColor(42);
241 pcharge1->SetXTitle("ADC units");
242 pcharge1->Draw();
243
244 pad13->cd();
245 xresid1->SetFillColor(42);
246 xresid1->Draw();
247
248 pad14->cd();
249 yresid1->SetFillColor(42);
250 yresid1->Draw();
251
252 TCanvas *c2 = new TCanvas("c2","Cluster Size",400,10,600,700);
253 pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99);
254 pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99);
255 pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49);
256 pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49);
257 pad21->SetFillColor(11);
258 pad22->SetFillColor(11);
259 pad23->SetFillColor(11);
260 pad24->SetFillColor(11);
261 pad21->Draw();
262 pad22->Draw();
263 pad23->Draw();
264 pad24->Draw();
265
266 pad21->cd();
267 npads1->SetFillColor(42);
268 npads1->SetXTitle("Cluster Size");
269 npads1->Draw();
270
271 pad23->cd();
272 xresys1->SetXTitle("x on pad");
273 xresys1->SetYTitle("x-xcog");
274 xresys1->Draw();
275
276 pad24->cd();
277 yresys1->SetXTitle("y on pad");
278 yresys1->SetYTitle("y-ycog");
279 yresys1->Draw();
280
281 TCanvas *c3 = new TCanvas("c3","Charge and Residuals",400,10,600,700);
282 pad31 = new TPad("pad31"," ",0.01,0.51,0.49,0.99);
283 pad32 = new TPad("pad32"," ",0.51,0.51,0.99,0.99);
284 pad33 = new TPad("pad33"," ",0.01,0.01,0.49,0.49);
285 pad34 = new TPad("pad34"," ",0.51,0.01,0.99,0.49);
286 pad31->SetFillColor(11);
287 pad32->SetFillColor(11);
288 pad33->SetFillColor(11);
289 pad34->SetFillColor(11);
290 pad31->Draw();
291 pad32->Draw();
292 pad33->Draw();
293 pad34->Draw();
294
295 pad31->cd();
296 ccharge2->SetFillColor(42);
297 ccharge2->SetXTitle("ADC units");
298 ccharge2->Draw();
299
300 pad32->cd();
301 pcharge2->SetFillColor(42);
302 pcharge2->SetXTitle("ADC units");
303 pcharge2->Draw();
304
305 pad33->cd();
306 xresid2->SetFillColor(42);
307 xresid2->Draw();
308
309 pad34->cd();
310 yresid2->SetFillColor(42);
311 yresid2->Draw();
312
313 TCanvas *c4 = new TCanvas("c4","Cluster Size",400,10,600,700);
314 pad41 = new TPad("pad41"," ",0.01,0.51,0.49,0.99);
315 pad42 = new TPad("pad42"," ",0.51,0.51,0.99,0.99);
316 pad43 = new TPad("pad43"," ",0.01,0.01,0.49,0.49);
317 pad44 = new TPad("pad44"," ",0.51,0.01,0.99,0.49);
318 pad41->SetFillColor(11);
319 pad42->SetFillColor(11);
320 pad43->SetFillColor(11);
321 pad44->SetFillColor(11);
322 pad41->Draw();
323 pad42->Draw();
324 pad43->Draw();
325 pad44->Draw();
326
327 pad41->cd();
328 npads2->SetFillColor(42);
329 npads2->SetXTitle("Cluster Size");
330 npads2->Draw();
331
332 pad43->cd();
333 xresys2->SetXTitle("x on pad");
334 xresys2->SetYTitle("x-xcog");
335 xresys2->Draw();
336
337 pad44->cd();
338 yresys2->SetXTitle("y on pad");
339 yresys2->SetYTitle("y-ycog");
340 yresys2->Draw();
341
342}