- Welding section on absorber side (LHCVC2C_001)
[u/mrichter/AliRoot.git] / MUON / MUONdoubles.C
CommitLineData
a9e2aefa 1#include "iostream.h"
2
3void MUONdoubles (Int_t evNumber1=0,Int_t evNumber2=0)
4{
5 {
6 hprof = new TProfile("hprof","Profile dmin vs r",24,0,120,0,50);
7 }
8 hdist1 = new TH1F("hdist1","distance",100,0,10);
9 hdist2 = new TH1F("hdist2","distance",100,0,10);
10 Float_t a1=3.7/TMath::Sqrt(4);
11 Float_t a2=26/TMath::Sqrt(4);
12 for (Int_t i=1; i<10000; i++) {
13 Float_t x1,x2,y1,y2;
14 x1 = (2*gRandom->Rndm(i)-1.)*a1;
15 x2 = (2*gRandom->Rndm(i)-1.)*a1;
16 y1 = (2*gRandom->Rndm(i)-1.)*a1;
17 y2 = (2*gRandom->Rndm(i)-1.)*a1;
18 Float_t d=TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
19 hdist1->Fill(d,1);
20 }
21 Float_t xh[100], yh[100];
22 for (Int_t k=0; k<1000; k++) {
23 for (Int_t i=0; i<100; i++) {
24 xh[i] = (2*gRandom->Rndm(i)-1.)*a2;
25 yh[i] = (2*gRandom->Rndm(i)-1.)*a2;
26 }
27 Float_t x1,y1;
28
29 for (Int_t i=0; i<10; i++) {
30 Float_t dmin=1000;
31 x1 = (2*gRandom->Rndm(i)-1.)*a2;
32 y1 = (2*gRandom->Rndm(i)-1.)*a2;
33 for (Int_t j=0; j<100; j++) {
34 Float_t d=TMath::Sqrt((x1-xh[j])*(x1-xh[j])+(y1-yh[j])*(y1-yh[j]));
35 if (d<dmin) dmin=d;
36 }
37 hdist2->Fill(dmin,1);
38 }
39 }
40
41
42
43/////////////////////////////////////////////////////////////////////////
44// This macro is a small example of a ROOT macro
45// illustrating how to read the output of GALICE
46// and do some analysis.
47//
48/////////////////////////////////////////////////////////////////////////
49
50// Dynamically link some shared libs
51
52 if (gClassTable->GetID("AliRun") < 0) {
53 gROOT->LoadMacro("loadlibs.C");
54 loadlibs();
55 }
56
57// Connect the Root Galice file containing Geometry, Kine and Hits
58
59 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
60 if (!file) file = new TFile("galice.root","UPDATE");
61
62// Get AliRun object from file or create it if not on file
63
64 if (!gAlice) {
65 gAlice = (AliRun*)file->Get("gAlice");
66 if (gAlice) printf("AliRun object found on file\n");
67 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
68 }
69//
70// Set reconstruction models
71//
72// Get pointers to Alice detectors and Digits containers
73 AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
74//
75// Book some histograms
76 TH1F *Hcentre[10], *Hall[10], *Hcfrac[10], *Htail[10], *Htfrac[10];
77 TH1F *Hdcut[10], *Hdall[10], *Hdfrac[10];
78
79 for (Int_t i=0; i<10; i++) {
80 Hcentre[i] = new TH1F("Hcentre","Hit Distribution",26,0,260);
81 Hcentre[i]->SetFillColor(0);
82 Hcentre[i]->SetXTitle("R (cm)");
83
84 Hall[i] = new TH1F("Hall","Hit Distribution",26,0,260);
85 Hall[i]->SetFillColor(0);
86 Hall[i]->SetXTitle("R (cm)");
87
88 Hcfrac[i] = new TH1F("Hcfrag","Hit Distribution",26,0,260);
89 Hcfrac[i]->SetFillColor(0);
90 Hcfrac[i]->SetXTitle("R (cm)");
91
92 Htail[i] = new TH1F("Htail","Hit Distribution",26,0,260);
93 Htail[i]->SetFillColor(0);
94 Htail[i]->SetXTitle("R (cm)");
95
96 Htfrac[i] = new TH1F("Htfrag","Hit Distribution",26,0,260);
97 Htfrac[i]->SetFillColor(0);
98 Htfrac[i]->SetXTitle("R (cm)");
99
100 Hdall[i] = new TH1F("Hdall","Hit Distribution",26,0,260);
101 Hdall[i]->SetFillColor(0);
102 Hdall[i]->SetXTitle("R (cm)");
103
104 Hdcut[i] = new TH1F("Hdcut","Hit Distribution",26,0,260);
105 Hdcut[i]->SetFillColor(0);
106 Hdcut[i]->SetXTitle("R (cm)");
107
108 Hdfrac[i] = new TH1F("Hdfrac","Hit Distribution",26,0,260);
109 Hdfrac[i]->SetFillColor(0);
110 Hdfrac[i]->SetXTitle("R (cm)");
111
112 }
113
114//
115// Loop over events
116//
117 Int_t Nh=0;
118 Int_t Nh1=0;
119 for (int nev=0; nev<= evNumber2; nev++) {
120 Int_t nparticles = gAlice->GetEvent(nev);
121 cout << "nev " << nev <<endl;
122 cout << "nparticles " << nparticles <<endl;
123 if (nev < evNumber1) continue;
124 if (nparticles <= 0) return;
125
126 TTree *TH = gAlice->TreeH();
127 Int_t ntracks = TH->GetEntries();
128 cout<<"ntracks "<<ntracks<<endl;
129
130 Int_t nbytes = 0;
131
132
133 TClonesArray *Particles = gAlice->Particles();
134 TTree *TD = gAlice->TreeD();
135 Int_t nent=TD->GetEntries();
136 printf("Found %d entries in the tree (must be one per cathode per event!)\n",nent);
137 Float_t x[2000];
138 Float_t y[2000];
139 Int_t nhit=0;
140
141 if (MUON) {
142//
143// Loop on chambers and on cathode planes
144//
145 TTree *TH = gAlice->TreeH();
146 Int_t ntracks = TH->GetEntries();
147
148//
149// Loop over Hits
150//
151 for (i=0; i<1000; i++) {
152 Float_t r =100.*gRandom->Rndm(i);
153 Float_t phi=2.*TMath::Pi()*gRandom->Rndm(i);
154 Float_t xr=r*TMath::Sin(phi);
155 Float_t yr=r*TMath::Cos(phi);
156 Hdall[0]->Fill(r,1.);
157 Float_t dmin=100000.;
158
159 if (i==0) {
160 for (Int_t track=0; track<ntracks; track++) {
161 gAlice->ResetHits();
162 Int_t nbytes += TH->GetEvent(track);
163 for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
164 mHit;
165 mHit=(AliMUONHit*)MUON->NextHit())
166 {
167 Int_t nch = mHit->fChamber; // chamber number
168 if (nch!=1) continue;
169
170 x[nhit] = mHit->fX; // x-pos of hit
171 y[nhit] = mHit->fY; // y-pos
172 nhit++;
173 } // hit loop
174 } // track loop
175 } else {
176 for (Int_t nh=0; nh<nhit; nh++) {
177 Float_t d=TMath::Sqrt((x[nh]-xr)*(x[nh]-xr)+(y[nh]-yr)*(y[nh]-yr));
178// printf ("\n r,d %f %f",r,d);
179 Float_t dx=TMath::Abs(x[nh]-xr);
180 Float_t dy=TMath::Abs(y[nh]-yr);
181 if (d < dmin) dmin=d;
182
183 if (dx<0.5 && dy < 0.75) Hdcut[0]->Fill(r,1.);
184 } // hit loop
185 } // i==0
186 if (r>20) hprof->Fill(r,dmin,1);
187 } // random loop
188 Int_t icat=0;
189 gAlice->ResetDigits();
190 gAlice->TreeD()->GetEvent(icat+1); // spurious +1 ...
191
192 for (Int_t ich=0;ich<1;ich++) {
193 AliMUONChamber iChamber= MUON->Chamber(ich);
194 TClonesArray *MUONdigits = MUON->DigitsAddress(ich);
195 if (MUONdigits == 0) continue;
196 //
197 // Get ready the current chamber stuff
198 //
199 AliMUONResponse* response = iChamber.GetResponseModel();
200 AliMUONSegmentation* seg = iChamber.GetSegmentationModel(icat+1);
201 HitMap = new AliMUONHitMapA1(seg, MUONdigits);
202 HitMap->FillHits();
203 Int_t nxmax=seg->Npx();
204 Int_t nymax=seg->Npy();
205//
206// generate random positions on the chamber
207//
208 for (Int_t i=0; i<4000; i++) {
209 Int_t ix = 2*Int_t((gRandom->Rndm(i)-0.5)*nxmax);
210 Int_t iy = 2*Int_t((gRandom->Rndm(i)-0.5)*nymax);
211
212// test for centre overlap
213//
214 Float_t xp,yp;
215 seg->GetPadCxy(ix,iy,xp,yp);
216 Float_t r=TMath::Sqrt(xp*xp+yp*yp);
217 if (TMath::Abs(xp)>3 && TMath::Abs(yp)>3)
218 Hall[ich]->Fill(r,1.);
219 if (HitMap->TestHit(ix,iy) != 0) {
220 Hcentre[ich]->Fill(r,1.);
221 }
222// test for neighbour overlap
223//
224 Int_t nn;
225 Int_t X[20], Y[20];
226 seg->Neighbours(ix,iy,&nn,X,Y);
227 Bool_t hit=false;
228
229 for (Int_t j=0; j<nn; j++) {
230 if (HitMap->TestHit(X[j],Y[j]) != 0) hit=true;
231 }
232 if (hit && HitMap->TestHit(ix,iy) == 0) Htail[ich]->Fill(r,1.);
233 } //random loop
234 delete HitMap;
235 } // chamber loop
236 } // event loop
237 } // if MUON
238
239 for (Int_t ich=0; ich<10; ich++) {
240 Hcfrac[ich]->Divide(Hcentre[ich],Hall[ich]);
241 Htfrac[ich]->Divide(Htail[ich],Hall[ich]);
242 Hdfrac[ich]->Divide(Hdcut[ich],Hdall[ich]);
243 }
244
245 TCanvas *c1 = new TCanvas("c1","Hit Densities",400,10,600,700);
246 pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
247 pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
248 pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
249 pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
250 pad11->SetFillColor(0);
251 pad12->SetFillColor(0);
252 pad13->SetFillColor(0);
253 pad14->SetFillColor(0);
254 pad11->Draw();
255 pad12->Draw();
256 pad13->Draw();
257 pad14->Draw();
258
259 pad11->cd();
260 Hcentre[0]->Draw();
261
262 pad12->cd();
263 Hall[0]->Draw();
264
265 pad13->cd();
266 Hcfrac[0]->Draw();
267
268 pad14->cd();
269 Htfrac[0]->Draw();
270
271 TCanvas *c2 = new TCanvas("c2","Hit Densities",400,10,600,700);
272 pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99);
273 pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99);
274 pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49);
275 pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49);
276 pad21->SetFillColor(0);
277 pad22->SetFillColor(0);
278 pad23->SetFillColor(0);
279 pad24->SetFillColor(0);
280 pad21->Draw();
281 pad22->Draw();
282 pad23->Draw();
283 pad24->Draw();
284
285 pad21->cd();
286 hdist1->Draw();
287
288 pad22->cd();
289 hdist2->Draw();
290
291 pad23->cd();
292 Hdfrac[0]->Draw();
293
294 pad24->cd();
295 hprof->Draw();
296 file->Close();
297}
298