]>
Commit | Line | Data |
---|---|---|
9c02cab7 | 1 | #include "PMD_dig_reco.h" |
2 | #include <iostream.h> // define cout stream | |
3 | #include <fstream.h> // defines ofstream class | |
4 | #include <stdlib.h> // defines exit() functions | |
5 | void PMD_dig_reco (Int_t evNumber=1) | |
6 | { | |
7 | ///////////////////////////////////////////////////////////////////////// | |
8 | // This macro is a small example of a ROOT macro | |
9 | // illustrating how to read the output of GALICE | |
10 | // and fill some histograms. | |
11 | // | |
12 | // | |
13 | // Created by : Baba Potukuchi.V.K.S. ( JU / CERN ) 17/10/2000. | |
14 | // Modified by : Baba Potukuchi.V.K.S. ( JU / CERN ) 6/02/2001. | |
15 | // e-mail : potukuchi.baba@cern.ch | |
16 | // | |
17 | ///////////////////////////////////////////////////////////////////////// | |
18 | // =========== Variables for Reconstruction============================== | |
19 | FILE *fp4; | |
20 | PMD_dig_reco* PPP = new PMD_dig_reco; | |
21 | // ====================================================================== | |
22 | Int_t padmult,cpvmult,padmult_gam,tot_trk,tot_trk_cpv; | |
23 | Int_t padmult_cut, cpvmult_cut; | |
24 | Int_t trk_Mult[100000]; | |
25 | Int_t trk_Mult_cpv[100000]; | |
26 | Float_t xPos,yPos,zPos,Edep,Vol1,Vol2,Vol2s,Vol3,Vol4,Vol5; | |
27 | Float_t xpad,ypad; | |
28 | Int_t iVol4,iVol5 , ixpad,iypad; | |
29 | Int_t nbytes = 0; | |
30 | // Float_t px,py,pz; | |
31 | Float_t px1,py1,pz1; | |
32 | Float_t mass; | |
33 | Float_t edep_tot; | |
34 | Int_t id_part; | |
35 | Int_t gam_in; | |
36 | Int_t jjj; | |
37 | Int_t pip_in,pim_in; | |
38 | TObjArray *points; | |
39 | long double r,eta,px,py,pz; | |
40 | Float_t phi,Eparticle,theta1,theta; | |
41 | FILE *fp1; | |
42 | FILE *fp2; | |
43 | ||
44 | fp1=fopen("kine_vertex.out","w"); | |
45 | fp2=fopen("kine_pmdhit.out","w"); | |
46 | fp4=fopen("det1.out","w"); | |
47 | ||
48 | // Dynamical link some shared libs | |
49 | if (gClassTable->GetID("AliRun") < 0) { | |
50 | gROOT->LoadMacro("loadlibs.C"); | |
51 | loadlibs(); | |
52 | } | |
53 | ||
54 | // Connect the Root Galice file containing Geometry, Kine and Hits | |
55 | TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root"); | |
56 | if (!file) file = new TFile("galice.root"); | |
57 | ||
58 | // Get AliRun object from file or create it if not on file | |
59 | // if (!gAlice) { | |
60 | gAlice = (AliRun*)file->Get("gAlice"); | |
61 | if (gAlice) printf("AliRun object found on file\n"); | |
62 | if (!gAlice) gAlice = new AliRun("gAlice","Alice test program"); | |
63 | // } | |
64 | // Just for test | |
65 | //Loop over the requested number of events | |
66 | for (Int_t iEvt = 0; iEvt < evNumber; iEvt++) { | |
67 | ||
68 | printf("evt no %d \n",iEvt); | |
69 | // End Just for test | |
70 | ||
71 | ||
72 | ||
73 | // Import the Kine and Hits Trees for the event evNumber in the file | |
74 | Int_t nparticles = gAlice->GetEvent(iEvt); | |
75 | ||
76 | printf ("N of particles %d\n",nparticles); | |
77 | ||
78 | if (nparticles <= 0) return; | |
79 | ||
80 | Int_t nevt = iEvt + 1; | |
81 | ||
82 | // ----********---- vertex track search BEGIN | |
83 | ||
84 | //Loop over the requested number of events | |
85 | // Import the Kine and Hits Trees for the event evNumber in the file | |
86 | ||
87 | TObjArray *Particles = gAlice->Particles(); | |
88 | Float_t phi; | |
89 | Float_t pl,cutmax,etamin,etamax,pmom,smin,smax; | |
90 | Float_t *pxyz; | |
91 | TObjArray *points; | |
92 | TParticle *particle; | |
93 | AliPoints *pm; | |
94 | // for (Int_t it = 0; it < 100000; it++)trk_Mult[it]=0.; | |
95 | ||
96 | for (Int_t in_track = 0; in_track < nparticles; in_track++) { | |
97 | // for (Int_t in_track = 0; in_track < 10; in_track++) { | |
98 | particle = (TParticle*)gAlice->Particle(in_track); | |
99 | if (!particle) continue; | |
100 | id_part = particle->GetPdgCode(); | |
101 | // mass = TDatabasePDG::Instance()->GetParticle(id_part)->Mass(); | |
102 | // printf("_____________--------%d %f\n",id_part,mass); | |
103 | ||
104 | px = particle->Px(); | |
105 | py = particle->Py(); | |
106 | pz = particle->Pz(); | |
107 | if(px !=0 || py!=0 || pz!=0) { | |
108 | if( (px!=0 || py!=0) && pz < 0 ){ | |
109 | // if(id_part == 22)gam_in++; | |
110 | // printf("%d\n",gam_in); | |
111 | r = TMath::Sqrt(px*px + py*py); | |
112 | theta = TMath::ATan2(r,TMath::Abs(pz)); | |
113 | phi = TMath::ATan2(py,px); | |
114 | phi = (phi*180.)/3.14; | |
115 | // if(theta < -1.5707961)printf( " THETA=%f\n",theta) ; | |
116 | if(theta > -1.5707961)eta = -TMath::Log(TMath::Tan(0.5*theta)); | |
117 | theta1=theta*180./3.14159; | |
118 | ||
119 | // if(TMath::Abs(eta) > 2.3 && TMath::Abs(eta) < 3.5){ | |
120 | if(id_part == 22)gam_in++; | |
121 | if(id_part == 211)pip_in++; | |
122 | if(id_part == -211)pim_in++; | |
123 | fprintf(fp1," %d %d %f %f %f %f %f %f \n",nevt,id_part,px,py,pz,eta,theta1,phi); | |
124 | // } | |
125 | // printf("%d %d %d\n",gam_in,pip_in,pim_in); | |
126 | // printf("%f %f %f %f %d\n",ptg,pxg,pyg,pzg,id_part); | |
127 | } | |
128 | } | |
129 | } | |
130 | ||
131 | // 1999 file-- Get pointers to Alice detectors and Hits containers | |
132 | //------********----- vertex track ENDs | |
133 | ||
134 | Float_t x,y,z,mass,e; | |
135 | Int_t nbytes = 0; | |
136 | Int_t j,hit,ipart; | |
137 | Int_t npmd; | |
138 | Int_t mbytes = 0; | |
139 | Int_t nhits,gam_in; | |
140 | Int_t ipmd,id_part; | |
141 | Int_t ypad1; | |
142 | Int_t mpmd; | |
143 | const Int_t Npad2 = 72; | |
144 | ||
145 | Int_t dVol; | |
146 | // Int_t iEvt; | |
147 | ||
148 | Int_t sector,plane; | |
149 | TParticle *particle; | |
150 | AliPMDhit *pmdHit; | |
151 | ||
152 | // Get pointers to Alice detectors and Hits containers | |
153 | AliDetector *PMD = gAlice->GetDetector("PMD"); | |
154 | // TClonesArray *fParticles = gAlice->Particles(); | |
155 | if (PMD) TClonesArray *PMDhits = PMD->Hits(); | |
156 | if (PMD) TObjArray *points = PMD->Points(); | |
157 | ||
158 | TTree *TH = gAlice->TreeH(); | |
159 | // TTree *TK = gAlice->TreeK(); | |
160 | ||
161 | Int_t ntracks = TH->GetEntries(); | |
162 | printf("No.of Tracksin the TreeH %d \n", ntracks); | |
163 | ||
164 | // =========== Initializing the array ============================ | |
165 | for (Int_t iDet = 1; iDet <= 2 ; ++iDet) { | |
166 | printf(" Initializing for Detector %d \n",iDet) ; | |
167 | for (Int_t iBox = 1; iBox <= 30 ; ++iBox) { | |
168 | for (Int_t iCol = 1; iCol <= 72 ; ++iCol) { | |
169 | for (Int_t iRow = 1; iRow <= 72 ; ++iRow) { | |
170 | if(iDet==1)(*PPP).module1 [iBox] [iCol] [iRow] = 0.0 ; | |
171 | if(iDet==2)(*PPP).module2 [iBox] [iCol] [iRow] = 0.0 ; | |
172 | } | |
173 | } | |
174 | } | |
175 | } | |
176 | // Creating histograms . | |
177 | TH1F *edp = new TH1F("edp","Energy Deposited in each cell" ,100,100.,100000.); | |
178 | TH1F *Modnum = new TH1F("Modnum","Module number" ,100,1.,100.); | |
179 | TH2F *xycells = new TH2F("xycells","XY-of Cells in an Unit Modlule " ,100,0.,80.,100,0.,80.); | |
180 | TH2F *xyhits = new TH2F("xyhits","Hits of the cells " ,100,-150.,150.,100,-150.,150.); | |
181 | ||
182 | // ================================================================ | |
183 | ||
184 | // ======================================================= | |
185 | // Start loop on tracks in the hits containers | |
186 | for (Int_t track=0; track<ntracks;track++) { | |
187 | gAlice->ResetHits(); | |
188 | ||
189 | //changes introduced ****** | |
190 | ||
191 | nbytes += TH->GetEvent(track); | |
192 | // mbytes += TK->GetEvent(track); | |
193 | // printf("%d ******** %d \n" ,track,nbytes); | |
194 | // ======>Histogram PMD | |
195 | ||
196 | if (PMD) { | |
197 | npmd = PMDhits->GetEntriesFast(); | |
198 | //printf("npmd=%d, for track=%d %d \n" ,npmd,track,nbytes); | |
199 | for (ipmd = 0; ipmd < npmd; ipmd++) { | |
200 | pmdHit = (AliPMDhit*) PMDhits->UncheckedAt(ipmd); | |
201 | ipart = pmdHit->GetTrack(); | |
202 | ||
203 | // printf("ID part and hits on PMD --------%d %d \n" ,ipart,track); | |
204 | // get kinematics of the particles | |
205 | // particle = (TParticle*)fParticles->UncheckedAt(ipart); | |
206 | particle = gAlice->Particle(ipart); | |
207 | // if (!particle) continue; | |
208 | // | |
209 | id_part = particle->GetPdgCode(); | |
210 | mass = TDatabasePDG::Instance()->GetParticle(id_part)->Mass(); | |
211 | // printf("ID part and hit on PMD %d %d \n" ,id_part,ipmd); | |
212 | px = particle->Px(); | |
213 | py = particle->Py(); | |
214 | pz = particle->Pz(); | |
215 | Eparticle = particle->Energy(); | |
216 | ||
217 | r = TMath::Sqrt(px*px + py*py); | |
218 | theta = TMath::ATan2(r,TMath::Abs(pz)); | |
219 | //phi = TMath::ATan2(TMath::Abs(py),TMath::Abs(px)); | |
220 | phi = TMath::ATan2(py,px); | |
221 | phi = (phi*180.)/3.14159; | |
222 | eta = -TMath::Log(TMath::Tan(0.5*theta)); | |
223 | if (pz < 0) eta = -eta; | |
224 | theta=theta*180./3.14159; | |
225 | ||
226 | // if(TMath::Abs(eta) > 2.3 && TMath::Abs(eta) < 3.5){ | |
227 | // | |
228 | xPos = pmdHit->X(); | |
229 | yPos = pmdHit->Y(); | |
230 | zPos = pmdHit->Z(); | |
231 | Vol1 = pmdHit->fVolume[0]; | |
232 | Vol2 = pmdHit->fVolume[1]; | |
233 | Vol3 = pmdHit->fVolume[2]; | |
234 | Vol4 = pmdHit->fVolume[3]; | |
235 | Vol5 = pmdHit->fVolume[4]; | |
236 | ||
237 | Vol2s = Vol2-1; | |
238 | ypad = Vol2s/Npad2 + 1; | |
239 | ypad1 = ypad-1; | |
240 | xpad = Vol2-ypad1*Npad2; | |
241 | Edep = pmdHit->fEnergy; | |
242 | ixpad = xpad+0.2 ; | |
243 | iypad = ypad+0.2 ; | |
244 | iVol4 = Vol4+0.2 ; | |
245 | iVol5 = Vol5+0.2 ; | |
246 | // module [Vol4] [Vol5] [xpad] [ypad] = Edep ; | |
247 | if(Edep>0)printf(" VOL4 %d %d %d %d %f \n",iVol4,iVol5,ixpad,iypad,Edep) ; | |
248 | edp ->Fill(Edep) ; | |
249 | Modnum ->Fill(iVol5) ; | |
250 | xycells->Fill(ixpad,iypad); | |
251 | xyhits->Fill(xPos,yPos); | |
252 | ||
253 | if(iVol4==1)(*PPP).module1 [iVol5] [ixpad] [iypad] = Edep ; | |
254 | if(iVol4==2)(*PPP).module2 [iVol5] [ixpad] [iypad] = Edep ; | |
255 | Float_t XXP = (*PPP).module2 [iVol5] [ixpad] [iypad] ; | |
256 | if(Edep>0)printf(" VOL4 %d %d %d %d %f \n",Vol4,Vol5,xpad,ypad,XXP) ; | |
257 | // Write to the kine.out file | |
258 | fprintf(fp2," %d %d %d %f %f %f %d %d %d %d %f %f %f %f %f \n", | |
259 | nevt,ipart,id_part,xPos,yPos,zPos,xpad,ypad,Vol5,Vol4,Edep,px,py,pz,Eparticle); | |
260 | ||
261 | // Closing Eta loop of the incident particles. | |
262 | // } | |
263 | ||
264 | } | |
265 | } | |
266 | } | |
267 | } | |
268 | ||
269 | pmdreco_(PPP); | |
270 | ||
271 | // TCanvas *c1 = new TCanvas("c1","PMD Simulated digits(Edep) ",400,10,600,700); | |
272 | // c1->Divide(1,2); | |
273 | // c1->cd(1); | |
274 | // gPad->SetFillColor(33); | |
275 | // edp->SetFillColor(42); | |
276 | // edp->Draw(); | |
277 | // c1->cd(2); | |
278 | // gPad->SetFillColor(33); | |
279 | // Modnum->SetFillColor(41); | |
280 | // Modnum->Draw(); | |
281 | // c1->Print("baba.ps"); | |
282 | // TCanvas *c2 = new TCanvas("c2","PMD Simulated digits(Edep) "); | |
283 | // gPad->SetFillColor(33); | |
284 | // xycells->SetFillColor(42); | |
285 | // xycells->Draw(); | |
286 | // c2->Print("baba1.ps"); | |
287 | // TCanvas *c3 = new TCanvas("c3","PMD Simulated digits(Edep) "); | |
288 | // gPad->SetFillColor(33); | |
289 | // xyhits->SetFillColor(42); | |
290 | // xyhits->Draw(); | |
291 | // c3->Print("baba2.ps"); | |
292 | } | |
293 | ||
294 | // =========== RECONSTRUCTION ====================================== | |
295 | void pmdreco_(PMD_dig_reco* PPP) { | |
296 | Float_t sqrth; | |
297 | sqrth = sqrt(3.) / 2.; | |
298 | ||
299 | printf(" I am in MAIN calling RanMAR & Intpts_ \n") ; | |
300 | rmarin_(PPP); | |
301 | intpts_(PPP); | |
302 | // Generate COORDINATES . | |
303 | for (Int_t i__ = 1; i__ <= 72; ++i__) { | |
304 | for (Int_t j = 1; j <= 72; ++j) { | |
305 | (*PPP).coord [1] [i__] [j] = i__ + j / 2.; | |
306 | (*PPP).coord [2] [i__] [j] = j * sqrth; | |
307 | } | |
308 | } | |
309 | //------------------------- | |
310 | for (Int_t iDet = 1; iDet <= 2 ; ++iDet) { | |
311 | printf("idet= %d \n",iDet) ; | |
312 | for (Int_t iBox = 1; iBox <= 30 ; ++iBox) { | |
313 | for (Int_t iCol = 1; iCol <= 72 ; ++iCol) { | |
314 | for (Int_t iRow = 1; iRow <= 72 ; ++iRow) { | |
315 | Float_t PVAL = (*PPP).module1 [iBox] [iCol] [iRow]; | |
316 | // printf("Det1 = %f %f %f %f %f \n",iDet, iBox, iCol, iRow,PVAL ) ; | |
317 | if(iDet==1)(*PPP).det1 [iCol] [iRow] = (*PPP).module1 [iBox] [iCol] [iRow] ; | |
318 | if(iDet==2)(*PPP).det2 [iCol] [iRow] = (*PPP).module2 [iBox] [iCol] [iRow] ; | |
319 | if(PVAL>0.0)printf("Det1 = %d %d %d %d %f \n",iDet, iBox, iCol, iRow,PVAL ) ; | |
320 | // fprintf(fp4," Det1 = %d %d %d %d %f \n",iDet, iBox, iCol, iRow,VAL ) ; | |
321 | ||
322 | } | |
323 | } | |
324 | // Ending Reading pads in the Box . | |
325 | } | |
326 | // Ending Reading pads for 1-CPV,2-PMD. | |
327 | } | |
328 | /* analyse detector 1 */ | |
329 | /* CPV clustering. Ordering data */ | |
330 | printf(" I am in MAIN Doing Clustering for CPV. \n") ; | |
331 | printf (" I am in MAIN calling ==order== for CPV(det=1) \n"); | |
332 | order_(1,PPP); | |
333 | printf(" I am in MAIN after order for CPV(det=1) \n") ; | |
334 | printf(" I am in MAIN calling ==crclust== \n") ; | |
335 | crclust_(1, PPP) ; | |
336 | printf(" I am in MAIN after ==crclust== & now call prints \n") ; | |
337 | prints_(1,PPP) ; | |
338 | printf(" I am in MAIN Doing Clustering for PMD. \n") ; | |
339 | printf (" I am in MAIN calling ==order== for PMD(det=2) \n"); | |
340 | order_(2,PPP); | |
341 | printf(" I am in MAIN after order for PMD(det=2) \n") ; | |
342 | printf(" I am in MAIN calling ==crclust== \n") ; | |
343 | crclust_(2, PPP) ; | |
344 | printf(" I am in MAIN after ==crclust== & now call prints \n") ; | |
345 | prints_(2,PPP) ; | |
346 | (*PPP).incr = 0; | |
347 | for (Int_t i__ = 1; i__ <= 72; ++i__) { | |
348 | for (j = 1; j <= 72; ++j) { | |
349 | if ((*PPP).infocl[1] [i__] [j] == 1) { | |
350 | ++(*PPP).incr; | |
351 | (*PPP).inford[1] [(*PPP).incr] = (*PPP).infocl[2] [i__] [j] ; | |
352 | (*PPP).inford[2] [(*PPP).incr] = i__; | |
353 | (*PPP).inford[3] [(*PPP).incr] = j; | |
354 | } | |
355 | } | |
356 | } | |
357 | for (i__ = 1; i__ <= 72; ++i__) { | |
358 | for (j = 1; j <= 72; ++j) { | |
359 | if ((*PPP).infocl[1] [i__] [j] == 2) { | |
360 | ++(*PPP).incr; | |
361 | (*PPP).inford[1] [(*PPP).incr] = (*PPP).infocl[2] [i__] [j]; | |
362 | (*PPP).inford[2] [(*PPP).incr] = i__; | |
363 | (*PPP).inford[3] [(*PPP).incr] = j; | |
364 | } | |
365 | } | |
366 | } | |
367 | arrange_(PPP); | |
368 | /* better clustering ( Gaussian fit ) */ | |
369 | refclust_(PPP); | |
370 | return ; | |
371 | } /* LOOP pmdreco */ | |
372 | // | |
373 | /* >>>>>>>>>>>> */void rmarin_(PMD_dig_reco* PPP) | |
374 | { | |
375 | Int_t i__, j, k, l, m; | |
376 | Float_t s, t; | |
377 | /* Initializing routine for RANMAR, must be called before */ | |
378 | /* generating any pseudorandom numbers with RANMAR. The input */ | |
379 | /* values should be in the ranges: 0<=IJ<=31328 */ | |
380 | /* 0<=KL<=30081 */ | |
381 | /* To get the standard values in Marsaglia's paper, IJ=1802, KL=9373 */ | |
382 | Int_t ij = 8775; | |
383 | Int_t kl = 5926; | |
384 | if (ij <= 0 || ij >= 31328) { | |
385 | goto L900; | |
386 | } | |
387 | if (kl <= 0 || kl >= 30081) { | |
388 | goto L900; | |
389 | } | |
390 | printf("I am in RanMAR \n"); | |
391 | i__ = ij / 177 % 177 + 2; | |
392 | j = ij % 177 + 2; | |
393 | k = kl / 169 % 178 + 1; | |
394 | l = kl % 169; | |
395 | printf(" ranmar initialized: %d %d %d %d %d %d \n",ij,kl,i__,j,k,l); | |
396 | for (Int_t ii = 1; ii <= 97; ++ii) { | |
397 | s = 0.; | |
398 | t = .5; | |
399 | for (Int_t jj = 1; jj <= 24; ++jj) { | |
400 | m = (((i__ * j) % 179) * k) % 179; | |
401 | i__ = j; | |
402 | j = k; | |
403 | k = m ; | |
404 | l = l * 53 + 1 % 169; | |
405 | if (l * m % 64 >= 32) { | |
406 | s += t; | |
407 | } | |
408 | t *= .5; | |
409 | } | |
410 | (*PPP).u[ii - 1] = s; | |
411 | } | |
412 | (*PPP).c__ = .021602869033813477; | |
413 | (*PPP).cd = .45623308420181274; | |
414 | (*PPP).cm = .99999982118606567; | |
415 | printf(" Exiting from rmarin %f %f %f \n",(*PPP).c__,(*PPP).cd,(*PPP).cm) ; | |
416 | return ; | |
417 | L900: | |
418 | printf( " +++++++ stop in rmarin: wrong params! "); | |
419 | } /* rmarin_ */ | |
420 | ||
421 | /*>>>>>>>>>>>>>>>> */void order_(Int_t jdet,PMD_dig_reco* PPP) | |
422 | { | |
423 | Int_t nmx,i__1, i__2, i__3; | |
424 | ||
425 | Float_t adum; | |
426 | Int_t idum, i__, j, k, l; | |
427 | Float_t d1[5185]; | |
428 | Float_t iord1[5185]; | |
429 | Int_t itest, i1, i2, id ; | |
430 | ||
431 | printf("I am in order for DET=%d \n",jdet); | |
432 | nmx = 5184; | |
433 | for (i1 = 1; i1 <= 72; ++i1) { | |
434 | for (i2 = 1; i2 <= 72; ++i2) { | |
435 | i__ = (i2 - 1) * 72 + i1; | |
436 | d1[i__] = 0.0 ; | |
437 | iord1 [i__] = 0; | |
438 | if(jdet==1) d1[i__] = (*PPP).det1 [i1] [i2] ; | |
439 | if(jdet==2) d1[i__] = (*PPP).det2 [i1] [i2] ; | |
440 | iord1 [i__] = i__; | |
441 | } | |
442 | } | |
443 | for (j = 2; j <= nmx; ++j) { | |
444 | itest = 0; | |
445 | adum = d1[j]; | |
446 | idum = iord1 [j]; | |
447 | i__2 = j - 1; | |
448 | for (k = 1; k <= i__2; ++k) { | |
449 | if (adum > d1[k ] && itest == 0) { | |
450 | itest = 1 ; | |
451 | i__3 = k; | |
452 | for (l = j - 1; l >= i__3; --l) { | |
453 | d1[l+1] = d1[l]; | |
454 | iord1 [l+1] = iord1 [l]; | |
455 | } | |
456 | d1[k] = adum; | |
457 | iord1 [k] = idum; | |
458 | } | |
459 | } | |
460 | } | |
461 | i__1 = nmx; | |
462 | for (i__ = 1; i__ <= i__1; ++i__) { | |
463 | id = iord1 [i__]; | |
464 | i2 = id / 72 + 1; | |
465 | i1 = id - (i2 - 1) * 72; | |
466 | (*PPP).iorder[1] [i__] = i1; | |
467 | (*PPP).iorder[2] [i__] = i2; | |
468 | ||
469 | printf("I am in order for DET=%d %d %d %d %d %d \n",jdet,i__,i1,i2,(*PPP).iorder[1] [i__],(*PPP).iorder[2] [i__]); | |
470 | } | |
471 | } /* order_ */ | |
472 | ||
473 | /* crclust_ */ void crclust_(Int_t ndet, PMD_dig_reco* PPP) | |
474 | { | |
475 | Float_t d__ [73] [73] ; | |
476 | ||
477 | Int_t neibx[7] = { 0,1,0,-1,-1,0,1 }; | |
478 | Int_t neiby[7] = { 0,0,1,1,0,-1,-1 }; | |
479 | ||
480 | Int_t i__1; | |
481 | Float_t sqrth ; | |
482 | ||
483 | Int_t icld, idum, i__, j, k; | |
484 | Int_t nb, id1, id2, jd1, jd2, icle, icn, itr, nmx ; | |
485 | ||
486 | printf(" neibx %d %d %d %d %d %d \n",neibx[1],neibx[2],neibx[3],neibx[4],neibx[5],neibx[6]) ; | |
487 | printf(" neiby %d %d %d %d %d %d \n",neiby[1],neiby[2],neiby[3],neiby[4],neiby[5],neiby[6]) ; | |
488 | ||
489 | for (j = 1; j <= 72; ++j) { | |
490 | for (k = 1; k <= 72; ++k) { | |
491 | d__ [j] [k] = 0.0 ; | |
492 | if(ndet==1) d__ [j] [k] = (*PPP).det1 [j] [k] ; | |
493 | if(ndet==2) d__ [j] [k] = (*PPP).det2 [j] [k] ; | |
494 | ||
495 | } | |
496 | } | |
497 | nmx = 5148; | |
498 | sqrth=TMath::Sqrt(3.)/2. ; | |
499 | /* initialise infocl */ | |
500 | for (i__ = 1; i__ <= 2; ++i__) { | |
501 | for (j = 1; j <= 72; ++j) { | |
502 | for (k = 1; k <= 72; ++k) { | |
503 | (*PPP).infocl[i__] [j] [k] = 0; | |
504 | } | |
505 | } | |
506 | } | |
507 | /* compute number of pads with nonzero count ( nmx1 ) */ | |
508 | i__1 = nmx; | |
509 | for (i__ = 1; i__ <= i__1; ++i__) { | |
510 | id1 = (*PPP).iorder [1] [i__]; | |
511 | id2 = (*PPP).iorder [2] [i__]; | |
512 | if (d__[id1] [id2] > 0.0000000000001) { | |
513 | (*PPP).nmx1 = i__; | |
514 | printf(" nmx1 -> %d %d %f %d \n",id1,id2, d__ [id1] [id2],(*PPP).nmx1 ) ; | |
515 | } | |
516 | } | |
517 | printf(" nmx1= %d \n",(*PPP).nmx1) ; | |
518 | /* crude clustering. start with pad with largest count ( cluster 1 ) | |
519 | */ | |
520 | icle = 1; | |
521 | idum = 0; | |
522 | id1 = (*PPP).iorder [1] [1]; | |
523 | id2 = (*PPP).iorder [2] [1]; | |
524 | (*PPP).infocl[1] [id1] [id2] = 1; | |
525 | (*PPP).infocl[2] [id1] [id2] = icle; | |
526 | printf(" infocl %d %d %d %d %d \n",id1,id2,icle,(*PPP).infocl[1] [id1] [id2],(*PPP).infocl[2] [id1] [id2] ) ; | |
527 | /* check if neighbouring count in nonzero. If yes, include in the */ | |
528 | /* cluster */ | |
529 | /* infocl(1,i1,i2) =1 for primary ( peak ) =2 for secondary ( */ | |
530 | /* neighbour ) */ | |
531 | /* infocl(2,i1,i2) stores cluster no. */ | |
532 | for (nb = 1; nb <= 6; ++nb) { | |
533 | jd1 = id1 + neibx[nb]; | |
534 | jd2 = id2 + neiby[nb]; | |
535 | printf(" infocl %d %d %d %d %d %f %d %d\n",nb,jd1,id1,jd2,id2,d__[jd1] [jd2],icle,(*PPP).infocl[1] [jd1] [jd2] ) ; | |
536 | if (jd1 > 0 && jd1 <= 72 && jd2 > 0 && jd2 <= 72) { | |
537 | if (d__[jd1] [jd2] > 0.0) { | |
538 | (*PPP).infocl[1] [jd1] [jd2] = 2; | |
539 | (*PPP).infocl[2] [jd1] [jd2] = icle; | |
540 | } | |
541 | } | |
542 | printf(" infocl %d %d %d %d %d %f %d %d\n",nb,jd1,id1,jd2,id2,d__[jd1] [jd2],icle,(*PPP).infocl[1] [jd1] [jd2] ) ; | |
543 | } | |
544 | /* start checking with next largest count and so on */ | |
545 | i__1 = (*PPP).nmx1; | |
546 | for (i__ = 2; i__ <= i__1; ++i__) { | |
547 | itr = 0; | |
548 | id1 = (*PPP).iorder[1] [i__]; | |
549 | id2 = (*PPP).iorder[2] [i__]; | |
550 | printf("iorder %d %d %d %d %d %d %d\n",i__,id1,id2,(*PPP).infocl [1] [id1] [id2],jd1,jd2,(*PPP).infocl [1] [jd1] [jd2] ) ; | |
551 | /* nonzero infocl(1,id1,id2) means the pad is included in earlier | |
552 | */ | |
553 | /* cluster ( itr = 1 ) Otherwise itr = 0. */ | |
554 | if ((*PPP).infocl [1] [id1] [id2] > 0) { | |
555 | itr = 1; | |
556 | } | |
557 | if (itr == 0) { | |
558 | for (nb = 1; nb <= 6; ++nb) { | |
559 | jd1 = id1 + neibx[nb ]; | |
560 | jd2 = id2 + neiby[nb ]; | |
561 | /* check if a neighbouring pad belongs to previous cluster | |
562 | . If yes, */ | |
563 | /* icn = 1. Otherwise icn = 0. */ | |
564 | icn = 0; | |
565 | if (jd1 > 0 && jd1 <= 72 && jd2 > 0 && jd2 <= 72) { | |
566 | if ((*PPP).infocl [1] [jd1] [jd2] != 0) { | |
567 | icn = 1; | |
568 | icld = (*PPP).infocl [2] [jd1] [jd2] ; | |
569 | } | |
570 | } | |
571 | } | |
572 | /* itr = 0 and icn = 0 means new cluster. Store the cluster in | |
573 | fo. in */ | |
574 | /* infocl. */ | |
575 | if (TMath::Abs(icn) <= 0.1) { | |
576 | ++icle; | |
577 | (*PPP).infocl[1] [id1] [id2] = 1; | |
578 | (*PPP).infocl[2] [id1] [id2] = icle; | |
579 | printf("ID1 icn %d %d %d %d %d\n",icn,icle,id1,id2,(*PPP).infocl[2] [id1] [id2]); | |
580 | for (nb = 1; nb <= 6; ++nb) { | |
581 | /* neighbours of the peak put with the peak. */ | |
582 | jd1 = id1 + neibx[nb ]; | |
583 | jd2 = id2 + neiby[nb ]; | |
584 | if (jd1 > 0 && jd1 <= 72 && jd2 > 0 && jd2 <= 72) { | |
585 | if (d__[jd1] [jd2] != 0 && (*PPP).infocl [1] [jd1] [jd2] == 0) { | |
586 | (*PPP).infocl[1] [jd1] [jd2] = 2; | |
587 | (*PPP).infocl[2] [jd1] [jd2] = icle; | |
588 | printf("JD1 icn %d %d %d %d %d\n",icn,icle,jd1,jd2,(*PPP).infocl[2] [jd1] [jd2]); | |
589 | } | |
590 | } | |
591 | } | |
592 | } else { | |
593 | /* itr = 0 and icn = 1 means the pad belongs to cluster ic | |
594 | ld */ | |
595 | (*PPP).infocl[1] [id1] [id2] = 2; | |
596 | (*PPP).infocl[2] [id1] [id2] = icld; | |
597 | for (nb = 1; nb <= 6; ++nb) { | |
598 | /* and also the neighbours. */ | |
599 | jd1 = id1 + neibx[nb ]; | |
600 | jd2 = id2 + neiby[nb ]; | |
601 | if (jd1 > 0 && jd1 <= 72 && jd2 > 0 && jd2 <= 72) { | |
602 | if (d__[jd1] [jd2] != 0. && (*PPP).infocl[1] [jd1] [jd2] == 0) { | |
603 | (*PPP).infocl[1] [jd1] [jd2] = 2; | |
604 | (*PPP).infocl[2] [jd1] [jd2] = icld; | |
605 | } | |
606 | } | |
607 | } | |
608 | } else { | |
609 | /* itr = 0, so not a new cluster */ | |
610 | icld = (*PPP).infocl[2] [id1] [id2] ; | |
611 | for (nb = 1; nb <= 6; ++nb) { | |
612 | /* neighbours put with the previous cluster. */ | |
613 | jd1 = id1 + neibx[nb ]; | |
614 | jd2 = id2 + neiby[nb ]; | |
615 | if (jd1 > 0 && jd1 <= 72 && jd2 > 0 && jd2 <= 72) { | |
616 | if (d__[jd1] [jd2] != 0. && (*PPP).infocl[1] [jd1] [jd2] == 0) { | |
617 | (*PPP).infocl[1] [jd1] [jd2] = 2; | |
618 | (*PPP).infocl[2] [jd1] [jd2] = icld; | |
619 | } | |
620 | } | |
621 | } | |
622 | } | |
623 | } | |
624 | } | |
625 | printf(" ndet=%d,nmx1= %d , icle= %d icld= %d \n",ndet,(*PPP).nmx1,icle,icld) ; | |
626 | ||
627 | /* L2: */ | |
628 | /* L3: */ | |
629 | return ; | |
630 | } /* crclust_ */ | |
631 | ||
632 | /* Subroutine */ void intpts_(PMD_dig_reco* PPP) | |
633 | { | |
634 | Int_t i__1, i__2, i__3; | |
635 | ||
636 | Int_t ixpt; | |
637 | Float_t x, y, dx, dy; | |
638 | Int_t ix, iy; | |
639 | Float_t yl; | |
640 | Int_t idx, idy, ipt; | |
641 | Float_t wtx, wty; | |
642 | ||
643 | ixpt = 2; | |
644 | idx = -1; | |
645 | ipt = 0; | |
646 | dx = .25; | |
647 | i__1 = ixpt; | |
648 | for (ix = -ixpt; ix <= i__1; ++ix) { | |
649 | x = ix * dx; | |
650 | idx = -idx; | |
651 | i__2 = ix / ixpt ; | |
652 | wtx = (3. - idx - (TMath::Abs(i__2) ))/ 3. * dx; | |
653 | yl = 1. / TMath::Sqrt(3.) * (1. - TMath::Abs(x)); | |
654 | dy = yl * 2. / 4.; | |
655 | idy = -1; | |
656 | i__2 = ixpt; | |
657 | for (iy = -ixpt; iy <= i__2; ++iy) { | |
658 | idy = -idy; | |
659 | y = iy * dy; | |
660 | i__3 = iy / ixpt ; | |
661 | wty = (3. - idy - (TMath::Abs(i__3))) / 3. * dy; | |
662 | ++ipt; | |
663 | (*PPP).pts[1] [ipt] = x; | |
664 | (*PPP).pts[2] [ipt] = y; | |
665 | (*PPP).wts[ipt] = wtx * wty; | |
666 | } | |
667 | } | |
668 | printf(" Exiting from intpts_ %f %f \n", ipt,(*PPP).wts[ipt] ); | |
669 | return ; | |
670 | } /* intpts_ */ | |
671 | ||
672 | /* prints_ */ void prints_(Int_t ndet, PMD_dig_reco* PPP) | |
673 | { | |
674 | Int_t i__, j; | |
675 | ||
676 | if (ndet == 1) { | |
677 | FILE *fp1; | |
678 | FILE *fp2; | |
679 | fp1=fopen("CPV_prime.out","w"); | |
680 | fp2=fopen("CPV_secondary.out","w"); | |
681 | } | |
682 | if (ndet == 2) { | |
683 | FILE *fp1; | |
684 | FILE *fp2; | |
685 | fp1=fopen("PMD_prime.out","w"); | |
686 | fp2=fopen("PMD_secondary.out","w"); | |
687 | } | |
688 | for (i__ = 1; i__ <= 72; ++i__) { | |
689 | for (j = 1; j <= 72; ++j) { | |
690 | if( (*PPP).infocl [1] [i__] [j] == 1) { | |
691 | fprintf(fp1," %f %f %d %d %d \n",(*PPP).coord [1] [i__] [j],(*PPP).coord [2] [i__] [j],i__,j,(*PPP).infocl [2] [i__] [j] ) ; | |
692 | } | |
693 | if( (*PPP).infocl [1] [i__] [j] == 2) { | |
694 | fprintf(fp2," %f %f %d %d %d \n",(*PPP).coord [1] [i__] [j],(*PPP).coord [2] [i__] [j],i__,j,(*PPP).infocl [2] [i__] [j] ) ; | |
695 | } | |
696 | } | |
697 | } | |
698 | return ; | |
699 | } /* prints_ */ | |
700 | ||
701 | /* Subroutine */void arrange_(PMD_dig_reco* PPP) | |
702 | { | |
703 | Int_t i__1, i__2, i__3; | |
704 | ||
705 | Int_t ihld1, ihld2, ihld3, i__, j, k, itest; | |
706 | ||
707 | i__1 = (*PPP).incr; | |
708 | for (i__ = 2; i__ <= i__1; ++i__) { | |
709 | itest = 0; | |
710 | ihld1 = (*PPP).inford [1] [i__]; | |
711 | ihld2 = (*PPP).inford [2] [i__]; | |
712 | ihld3 = (*PPP).inford [3] [i__]; | |
713 | i__2 = i__ - 1; | |
714 | for (j = 1; j <= i__2; ++j) { | |
715 | if (itest == 0 && (*PPP).inford [1] [j] > ihld1) { | |
716 | itest = 1; | |
717 | i__3 = j; | |
718 | for (k = i__ -1; k >= i__3; --k) { | |
719 | (*PPP).inford [1] [k+1] = (*PPP).inford [1] [k] ; | |
720 | (*PPP).inford [2] [k+1] = (*PPP).inford [2] [k] ; | |
721 | (*PPP).inford [3] [k+1] = (*PPP).inford [3] [k] ; | |
722 | } | |
723 | (*PPP).inford [1] [j] = ihld1 ; | |
724 | (*PPP).inford [2] [j] = ihld2 ; | |
725 | (*PPP).inford [3] [j] = ihld3 ; | |
726 | } | |
727 | } | |
728 | } | |
729 | printf(" Exiting from arrange \n") ; | |
730 | return ; | |
731 | } /* arrange_ */ | |
732 | ||
733 | /* refclust_ */void refclust_(PMD_dig_reco* PPP) | |
734 | { | |
735 | ||
736 | Int_t i__1, i__2; | |
737 | Float_t d__ [73] [73] ; | |
738 | Int_t npad[201], i__, j, k; | |
739 | Int_t id, id1, id2, iclb; | |
740 | ||
741 | /* using information from crude clustering, a refined algorithm for */ | |
742 | /* clustering tried. At the moment using Gaussian fit with a number */ | |
743 | /* of Gaussian functions. */ | |
744 | for (j = 1; j <= 72; ++j) { | |
745 | for (k = 1; k <= 72; ++k) { | |
746 | d__ [j] [k] = (*PPP).det2 [j] [k] ; | |
747 | } | |
748 | } | |
749 | ||
750 | (*PPP).clno = 0; | |
751 | /* The No.of pads in each cluster = 200 */ | |
752 | for (i__ = 1; i__ <= 200; ++i__) { | |
753 | npad[i__] = 0; | |
754 | } | |
755 | for (i__ = 1; i__ <= 5000; ++i__) { | |
756 | (*PPP).ncl[i__] = 0; | |
757 | } | |
758 | /* number of pads in each cluster determined. */ | |
759 | i__1 = (*PPP).incr; | |
760 | for (i__ = 1; i__ <= i__1; ++i__) { | |
761 | id = (*PPP).inford[1] [i__]; | |
762 | ++(*PPP).ncl[id]; | |
763 | Int_t nclid = (*PPP).ncl[id] ; | |
764 | printf("id,ncl= %d %d\n",id,(*PPP).ncl[id]) ; | |
765 | } | |
766 | iclb = 0; | |
767 | /* number of clusters determined. */ | |
768 | for (i__ = 1; i__ <= 5000; ++i__) { | |
769 | id = (*PPP).ncl[i__]; | |
770 | if (id != 0 ) { | |
771 | ++iclb; | |
772 | ++npad[id]; | |
773 | } | |
774 | } | |
775 | printf("iclb= %d \n",iclb) ; | |
776 | id = 0; | |
777 | /* define coordinates of each pad in a cluster and call gausspara */ | |
778 | /* to determine Gaussian parameters. ( should it be eliminated ? ) */ | |
779 | i__1 = iclb; | |
780 | for (i__ = 1; i__ <= i__1; ++i__) { | |
781 | i__2 = (*PPP).ncl[i__]; | |
782 | printf("(*PPP).ncl(i)= %d %d\n",i__,(*PPP).ncl[i__]) ; | |
783 | for (j = 1; j <= i__2; ++j) { | |
784 | ++id; | |
785 | id1 = (*PPP).inford[2] [id]; | |
786 | id2 = (*PPP).inford[3] [id]; | |
787 | (*PPP).x[j] = (*PPP).coord[1] [id1] [id2]; | |
788 | (*PPP).y[j] = (*PPP).coord[2] [id1] [id2]; | |
789 | (*PPP).z__[j] = d__[id1] [id2]; | |
790 | } | |
791 | printf(" I am in refclust calling ==gausspara_== \n") ; | |
792 | gausspara_(i__2,PPP ); | |
793 | printf(" I am in refclust after ==gausspara_== \n") ; | |
794 | } | |
795 | printf("clno= %d \n",(*PPP).clno); | |
796 | return ; | |
797 | } /* refclust_ */ | |
798 | ||
799 | /* gausspara_ */ void gausspara_(Int_t n,PMD_dig_reco* PPP) | |
800 | { | |
801 | Int_t i__1, i__2, i__3; | |
802 | Float_t d__1, d__2; | |
803 | ||
804 | Int_t ihld, iord[501], idum, i__, j, k; | |
805 | Int_t itest; | |
806 | Int_t ig; | |
807 | Float_t pi, xi, yi, xj, yj, rrr; | |
808 | Float_t sum; | |
809 | ||
810 | /* single Gaussian of the form exp(-((x-x_0)^2+(y-y_0)^2)/r_0^2) */ | |
811 | ||
812 | printf(" I am in ==gausspara_== %d \n" , n) ; | |
813 | pi = 3.141593; | |
814 | Int_t nn = n ; | |
815 | if (nn == 1) { | |
816 | ++(*PPP).clno; | |
817 | /* single pad cluster. Assumed fit to single Gaussian centered at the */ | |
818 | /* center of the pad with Gaussian parameter r0 about 0.5 since */ | |
819 | /* neighbouring pads don't fire. */ | |
820 | (*PPP).x0[1] = (*PPP).x[1]; | |
821 | (*PPP).y0[1] = (*PPP).y[1]; | |
822 | (*PPP).r0[1] = 0.5; | |
823 | /* Computing 2nd power */ | |
824 | d__1 = (*PPP).r0[1]; | |
825 | (*PPP).h0[1] = (*PPP).z__[1] / pi / (d__1 * d__1); | |
826 | printf("single pad ,clust no. = %d %d %f %f \n",n, (*PPP).clno, (*PPP).x0[1],(*PPP).y0[1] ) ; | |
827 | } else if (nn == 2) { | |
828 | ++((*PPP).clno); | |
829 | /* two pad cluster. Assumed fit to single Gaussian centered at the */ | |
830 | /* weighted average of the centers of two pads with r_0 about 0.5. */ | |
831 | /* Could be a pair of clusters but no way of differentiating the two. */ | |
832 | /* print *,'two pads ','clust no. ', clno */ | |
833 | printf("two pad clusters = %d \n",(*PPP).clno ) ; | |
834 | // printf("z__(1),(2) = %f %f \n",(*PPP).z__[1],(*PPP).z__[2] ) ; | |
835 | Int_t btest = ((*PPP).z__[1])+((*PPP).z__[2]) ; | |
836 | if ( btest != 0 ) { | |
837 | (*PPP).x0[1] = ((*PPP).z__[1] * (*PPP).x[1] + (*PPP).z__[2] * (*PPP).x[2]) / ((*PPP).z__[1] + (*PPP).z__[2]); | |
838 | (*PPP).y0[1] = ((*PPP).z__[1] * (*PPP).y[1] + (*PPP).z__[2] * (*PPP).y[2]) / ((*PPP).z__[1] + (*PPP).z__[2]); | |
839 | (*PPP).r0[1] = 0.5; | |
840 | } | |
841 | /* Computing 2nd power */ | |
842 | d__1 = (*PPP).r0[1] = 0.5; | |
843 | (*PPP).h0[1] = ((*PPP).z__[1] + (*PPP).z__[2]) / pi / (d__1 * d__1); | |
844 | // printf("Two pad clusters = %d %d %f %f \n",n, (*PPP).clno, (*PPP).x0[1],(*PPP).y0[1] ) ; | |
845 | } else { | |
846 | /* guessing possible no. of Gaussians required. */ | |
847 | /* assuming that if pads are separated by more than one unit, they */ | |
848 | /* are possible candidates for Gaussian centers. */ | |
849 | /* do this from larger values of energy ( z here ). */ | |
850 | i__1 = nn; | |
851 | for (i__ = 1; i__ <= i__1; ++i__) { | |
852 | iord[i__] = i__; | |
853 | } | |
854 | i__1 = nn; | |
855 | for (i__ = 2; i__ <= i__1; ++i__) { | |
856 | itest = 0; | |
857 | ihld = iord[i__]; | |
858 | i__2 = i__ - 1; | |
859 | for (j = 1; j <= i__2; ++j) { | |
860 | if (itest == 0 && (*PPP).z__[iord[j]] < (*PPP).z__[ihld]) { | |
861 | itest = 1; | |
862 | i__3 = j; | |
863 | for (k = i__ - 1; k >= i__3; --k) { | |
864 | iord[k+1] = iord[k]; | |
865 | } | |
866 | iord[j] = ihld; | |
867 | } | |
868 | } | |
869 | } | |
870 | /* Gaussian centers determined. */ | |
871 | ig = 1; | |
872 | (*PPP).icl[ig] = iord[ig]; | |
873 | i__1 = nn; | |
874 | for (i__ = 2; i__ <= i__1; ++i__) { | |
875 | xi = (*PPP).x[iord[i__]]; | |
876 | yi = (*PPP).y[iord[i__]]; | |
877 | idum = 0; | |
878 | i__2 = ig; | |
879 | for (j = 1; j <= i__2; ++j) { | |
880 | xj = (*PPP).x[(*PPP).icl[j]]; | |
881 | yj = (*PPP).y[(*PPP).icl[j]]; | |
882 | /* Computing 2nd power */ | |
883 | d__1 = xi - xj; | |
884 | /* Computing 2nd power */ | |
885 | d__2 = yi - yj; | |
886 | rrr = TMath::Sqrt(d__1 * d__1 + d__2 * d__2); | |
887 | if (1.01 < rrr) { | |
888 | ++idum; | |
889 | } | |
890 | } | |
891 | if (idum == ig) { | |
892 | ++ig; | |
893 | (*PPP).icl[ig] = iord[i__]; | |
894 | } | |
895 | } | |
896 | (*PPP).clno += ig; | |
897 | printf("No. of pads = %d , No. of Gaussians = %f \n ",ig,(*PPP).clno ) ; | |
898 | sum = 0.; | |
899 | i__1 = ig; | |
900 | for (i__ = 1; i__ <= i__1; ++i__) { | |
901 | (*PPP).x0[i__] = (*PPP).x[(*PPP).icl[i__]]; | |
902 | (*PPP).y0[i__] = (*PPP).y[(*PPP).icl[i__]]; | |
903 | (*PPP).r0[i__] = .5; | |
904 | sum += (*PPP).z__[(*PPP).icl[i__]]; | |
905 | } | |
906 | i__1 = ig; | |
907 | for (i__ = 1; i__ <= i__1; ++i__) { | |
908 | if ( sum != 0 ) { | |
909 | (*PPP).h0[i__] = (*PPP).z__[(*PPP).icl[i__]] / sum; | |
910 | } | |
911 | } | |
912 | printf(" I am in ==gausspara_== Calling gausfit_ \n") ; | |
913 | gausfit_(ig, nn, PPP ); | |
914 | printf(" I am in ==gausspara_ after gausfit_ \n") ; | |
915 | } | |
916 | return ; | |
917 | } /* gausspara_ */ | |
918 | ||
919 | /* gausfit_ */ void gausfit_(Int_t ng, Int_t nn, PMD_dig_reco* PPP ) | |
920 | { | |
921 | Int_t i__1; | |
922 | Float_t d__1; | |
923 | ||
924 | Int_t iran; | |
925 | Float_t aint; | |
926 | Int_t i__; | |
927 | Float_t s; | |
928 | Float_t pi, ss, st; | |
929 | ||
930 | pi = 3.141593; | |
931 | printf(" ng=%d , nn=%d \n",ng,nn); | |
932 | s = 0.; | |
933 | i__1 = nn; | |
934 | for (i__ = 1; i__ <= i__1; ++i__) { | |
935 | s +=(*PPP).z__[i__]; | |
936 | } | |
937 | /* print *,s */ | |
938 | printf(" s=%d \n",s); | |
939 | i__1 = ng; | |
940 | for (i__ = 1; i__ <= i__1; ++i__) { | |
941 | /* Computing 2nd power */ | |
942 | // printf("i__, X0,y0,r0 = %d %f %f %f \n",i__, (*PPP).x0[i__],(*PPP).y0[i__],(*PPP).r0[i__]) ; | |
943 | d__1 = (*PPP).r0[i__]; | |
944 | (*PPP).h0[i__] = (*PPP).h0[i__] * s / pi / (d__1 * d__1); | |
945 | } | |
946 | st = 0.; | |
947 | i__1 = nn; | |
948 | for (i__ = 1; i__ <= i__1; ++i__) { | |
949 | Int_t pbab = 1 ; | |
950 | aintg_(ng, i__, pbab, aint, PPP); | |
951 | /* Computing 2nd power */ | |
952 | d__1 = aint - (*PPP).z__[i__]; | |
953 | if( s > 0. ) st += d__1 * d__1 / s; | |
954 | } | |
955 | // printf(" st=%d ",st); | |
956 | for (iran = 1; iran <= 5000; ++iran) { | |
957 | Float_t rand; | |
958 | i__1 = ng; | |
959 | for (i__ = 1; i__ <= i__1; ++i__) { | |
960 | ranmar_(rand ,PPP) ; | |
961 | (*PPP).xx[i__] = (*PPP).x0[i__] + (rand - .5) * 1.; | |
962 | ranmar_(rand ,PPP) ; | |
963 | (*PPP).yy[i__] = (*PPP).y0[i__] + (rand - .5) * 1.; | |
964 | ranmar_(rand ,PPP) ; | |
965 | (*PPP).rr[i__] = (*PPP).r0[i__] * (((rand -.5) * .2 )+ 1.); | |
966 | /* Computing 2nd power */ | |
967 | d__1 = (*PPP).r0[i__] / (*PPP).rr[i__]; | |
968 | (*PPP).hh[i__] = (*PPP).h0[i__] * (d__1 * d__1); | |
969 | } | |
970 | ss = 0.; | |
971 | i__1 = nn; | |
972 | for (i__ = 1; i__ <= i__1; ++i__) { | |
973 | Int_t pbab = 2 ; | |
974 | aintg_(ng, i__, pbab, aint, PPP); | |
975 | /* Computing 2nd power */ | |
976 | d__1 = aint - (*PPP).z__[i__]; | |
977 | if(s>0.0)ss += d__1 * d__1 / s; | |
978 | } | |
979 | if (ss < st) { | |
980 | i__1 = ng; | |
981 | for (i__ = 1; i__ <= i__1; ++i__) { | |
982 | (*PPP).x0[i__] = (*PPP).xx[i__]; | |
983 | (*PPP).y0[i__] = (*PPP).yy[i__]; | |
984 | (*PPP).r0[i__] = (*PPP).rr[i__]; | |
985 | (*PPP).h0[i__] = (*PPP).hh[i__]; | |
986 | } | |
987 | st = ss; | |
988 | printf(" st=%d \n",st); | |
989 | } | |
990 | } | |
991 | /* print 1,((*PPP).x0(i),(*PPP).(*PPP).y0(i),(*PPP).r0(i),i=1,ng) */ | |
992 | /* print *,st */ | |
993 | i__1 = ng; | |
994 | for (i__ = 1; i__ <= i__1; ++i__) { | |
995 | // printf("x0(i),y0(i),r0(i) =%f %f %f \n",(*PPP).x0[i__1],(*PPP).y0[i__1],(*PPP).r0[i__1] ) ; | |
996 | } | |
997 | printf(" I am ending gausfit_ \n") ; | |
998 | return ; | |
999 | } /* gausfit_ */ | |
1000 | ||
1001 | /* aintg_ */void aintg_(Int_t ng, Int_t n, Int_t bab, Float_t aint, PMD_dig_reco* PPP ) | |
1002 | { | |
1003 | Int_t i__1; | |
1004 | Float_t d__1, d__2, d__3; | |
1005 | ||
1006 | Int_t i__, j; | |
1007 | Float_t xp, yp, fun; | |
1008 | ||
1009 | aint = 0.; | |
1010 | for (i__ = 1; i__ <= 25; ++i__) { | |
1011 | xp = (*PPP).pts[1] [i__] + (*PPP).x[n] ; | |
1012 | yp = (*PPP).pts[2] [i__] + (*PPP).y[n] ; | |
1013 | // printf("x[n],y[n],pts[1] [i__],[2] [i__], xp,yp =%f %f %f %f %f %f\n",(*PPP).x[n],(*PPP).y[n],(*PPP).pts[1] [i__],(*PPP).pts[2] [i__],xp,yp) ; | |
1014 | fun = 0.; | |
1015 | i__1 = ng; | |
1016 | for (j = 1; j <= i__1; ++j) { | |
1017 | //printf("J,D__1,D__2,D__3= %d %f %f %f %f\n",j,(*PPP).xx[j],xp,(*PPP).yy[j],(*PPP).rr[j]) ; | |
1018 | /* Computing 2nd power */ | |
1019 | d__1 = xp - (*PPP).xx[j]; | |
1020 | /* Computing 2nd power */ | |
1021 | d__2 = yp - (*PPP).yy[j]; | |
1022 | /* Computing 2nd power */ | |
1023 | if( bab == 1) d__3 = (*PPP).r0[j]; | |
1024 | if( bab == 2) d__3 = (*PPP).rr[j]; | |
1025 | // printf("J,D__1,D__2,D__3= %d %f %f %f \n",j,d__1,d__2,d__3) ; | |
1026 | fun += TMath::exp(-(d__1 * d__1 + d__2 * d__2) / (d__3 * d__3)) * (*PPP).hh[j]; | |
1027 | } | |
1028 | aint += fun * (*PPP).wts[i__]; | |
1029 | // printf("FUN,aint= %f %f \n",fun,aint) ; | |
1030 | } | |
1031 | return ; | |
1032 | } /* aintg_ */ | |
1033 | /* ranmar() >*/void ranmar_(Float_t random , PMD_dig_reco* PPP) | |
1034 | { | |
1035 | /* Initialized data */ | |
1036 | ||
1037 | Int_t i__ = 97; | |
1038 | Int_t j = 33; | |
1039 | Float_t uni; | |
1040 | ||
1041 | /* ===========================================C==========================C | |
1042 | */ | |
1043 | /* Universal random number generator proposed by Marsaglia and Zaman */ | |
1044 | /* in report FSU-SCRI-87-50 */ | |
1045 | uni = (*PPP).u[i__ - 1] - (*PPP).u[j - 1]; | |
1046 | if (uni < 0.) { | |
1047 | uni += 1.; | |
1048 | } | |
1049 | (*PPP).u[i__ ] = uni; | |
1050 | --i__; | |
1051 | if (i__ == 0) { i__ = 97; } | |
1052 | --j; | |
1053 | if (j == 0) { j = 97; } | |
1054 | (*PPP).c__ -= (*PPP).cd; | |
1055 | if ((*PPP).c__ < 0.) { (*PPP).c__ += (*PPP).cm; } | |
1056 | uni -= (*PPP).c__; | |
1057 | if (uni <= 0.) { uni += 1.; } | |
1058 | random = uni; | |
1059 | } /* ranmar_ */ | |
1060 |