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